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SUMMARY 

We  present  a  new  method  for  evaluating  the  inner  integral  of  the  impedance  matrix  element  in 
the  traditional  Rao-Wilton-Glisson  formulation  of  the  method  of  moments  for  perfect 
conductors.  In  this  method  we  replace  the  original  integrand  (modified  by  a  constant  phase 
factor)  by  its  Taylor  series  and  keep  enough  terms  to  guarantee  a  number  of  significant  digits  in 
the  integration  outcome.  We  develop  criteria  that  relate  the  number  of  Taylor  tenns  to  the 
number  of  required  significant  digits.  We  integrate  the  leading  Taylor  terms  analytically  and  the 
rest  through  iteration  formulas.  We  show  that  the  iteration  fonnulas  converge  for  all  observation 
points  within  a  sphere  with  a  radius  of  half-a-wavelength  and  center  the  triangle’s  centroid.  We 
compare  results  of  our  method  with  existing  ones  and  find  them  in  excellent  agreement.  Outside 
this  sphere,  we  employ  a  set  of  triangle  cubatures  of  increasing  size  together  with  a  convergence 
criterion  to  determine  the  integration  outcome  to  a  prescribed  number  of  significant  digits.  By 
designing  appropriate  numerical  experiments  using  a  set  of  25  triangles  and  about  10,000 
observation  points,  we  define  spherical  regions  of  space  where  a  cubature  of  minimum  size  will 
provide  a  desired  number  of  significant  digits.  The  approach  is  quite  general  but  we  demonstrate 
it  explicitly  by  using  seven  significant  digits  as  the  required  accuracy. 
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INTRODUCTION 

The  Method  of  Moments  (MoM)  is  one  of  the  principal  methods  for  electromagnetic 
simulations  in  the  frequency  domain.  In  this  method,  we  express  an  unknown  current  density  as  a 
linear  combination  of  known  functions  and  the  objective  of  the  MoM  is  to  detennine  the 
coefficients  in  this  expansion  by  minimizing  the  square  of  the  modulus  of  the  residual  error.  It 
ultimately  leads  to  a  system  of  linear  algebraic  equations,  the  solution  of  which  yields  the  values 
of  the  coefficients  of  the  current  density  expansion.  The  elements  of  the  matrix  (known  as  the 
impedance  matrix  (IM))  are  expressed  as  double  surface  integrals  over  flat,  triangular  regions  of 
space.  The  purpose  of  this  study  is  to  compute  the  inner  of  these  integrals  to  a  prescribed 
precision  when  the  object  under  consideration  is  a  perfect  conductor. 

Until  recently,  the  quantities  of  interest  in  electromagnetic  simulations  were  the  far-field 
electric  and  magnetic  fields.  From  these,  we  could  compute  the  radar  cross  section  of  a  target  or 
the  radiation  patterns  of  an  antenna.  This  allowed  for  a  considerable  degree  of  error  in  computing 
the  elements  of  the  IM  in  the  MoM  because  of  the  error  smoothing  effect  of  the  near-  to  far-field 
transfonnation  (integration).  Moreover,  due  to  computer  hardware  limitations,  the  size  of  the  IM 
was  small  enough  so  that  the  round-off  error  did  not  have  a  severe  effect  on  the  accuracy  of  the 
solution. 

In  recent  years,  the  MoM  is  being  applied  to  problems  where  near-field  infonnation  is 
required  (e.g.,  antennas  and  antenna  arrays).  This  requires  a  more  accurate  computation  of  the 
elements  of  the  IM  than  when  only  the  far  fields  are  of  interest.  Additionally,  advances  in 
computer  hardware  allow  us  to  address  problems  that  result  in  an  IM  system  with  millions  of 
unknowns.  Thus,  the  effect  of  the  round-off  error  becomes  more  pronounced.  Both  of  these 
reasons  lead  us  to  the  conclusion  that  the  more  accurate  the  representation  of  the  IM  is,  the  better 
the  quality  of  the  solution.  This  is  the  central  theme  of  the  present  study:  the  computation  of  the 
inner  surface  integral  of  an  IM  element  to  prescribed  precision  when  the  platform  of  interest  is 
perfectly  conducting. 

The  calculation  of  the  inner  surface  integral  comprises  two  parts.  In  the  first  part  we  present 
an  approximation  method  based  on  Taylor’s  Theorem  with  a  Remainder.  This  theorem  allows  us 
to  approximate  the  integrand  (that  we  do  not  know  how  to  integrate  analytically)  by  a 
polynomial  that  we  can  integrate  analytically.  The  degree  of  approximation,  i.e.,  the  number  of 
significant  digits  to  which  the  approximation  is  accurate  is  known  and,  hence,  so  is  that  of  the 
integration.  The  answer  to  this  evaluation  is  given  in  tenns  of  elementary  transcendental 
functions  for  which  there  exist  robust  computational  algorithms.  Where  these  algorithms  might 
fail,  we  have  replaced  them  by  approximations  that  yield  the  required  number  of  significant 
digits. 

The  drawback  in  using  this  method  to  compute  the  inner  integral  of  the  IM  element  is  that  the 
process  that  involves  the  analytical  evaluation  of  the  approximated  surface  integral  is  based  on 
two  iteration  formulas  that  do  not  converge  for  all  values  of  the  observation  point.  This  point  is 
the  integration  variable  of  the  outer  surface  integral  of  the  IM  element.  We  have  shown, 
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however,  that  the  iteration  formulas  do  converge  for  all  values  of  the  observation  point  within  a 
distance  of  half  a  wavelength  from  the  centroid  of  the  triangle.  Outside  this  region,  we  show  in 
the  second  part  of  this  report  that  we  can  use  existing  cubatures1  to  introduce  spherical  regions  in 
the  observation-point  space  where  a  cubature  of  a  minimum  size  will  produce  the  required 
number  of  significant  digits.  These  spherical  regions  have  their  center  at  the  triangle’s  centroid. 

In  conclusion,  we  have  developed  a  procedure  whereby  we  divide  the  entire  space  around  a 
triangle  into  spherical  regions  with  center  the  triangle’s  centroid.  In  the  innermost  region,  we  use 
approximation  theory  to  determine  the  value  of  the  inner  integral  of  the  IM  element  to  a 
prescribed  number  of  significant  digits.  In  this  region,  there  is  no  cubature  that  can  yield  the 
same  precision  in  a  shorter  time.  In  the  rest  of  the  spherical  regions,  we  use  the  limit  of  a  (finite) 
sequence  of  cubatures  to  determine  the  value  of  the  inner  integral  to  the  prescribed  number  of 
significant  digits.  Subsequently,  in  each  spherical  region,  we  use  the  smallest  cubature  that  will 
yield  this  accuracy.  We  have  found  that  three  spherical  regions  are  sufficient  to  cover  the  entire 
observation  space. 

We  know  of  no  other  method  that  calculates  the  inner  integral  of  the  IM  element  to  a 
prescribed  precision. 


1  By  cubature  we  mean  a  numerical  algorithm  that  evaluates  a  surface  integral.  We  reserve  the  term  quadrature  for 
an  algorithm  that  evaluates  a  line  integral. 
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PART  1.  ANALYTICAL  EVALUATION  OF  INTEGRALS 
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1.  STATEMENT  OF  PROBLEM  AND  APPROACH 

In  the  Rao-Wilton-Glisson  [JJ  formulation  of  the  method  of  moments  (MoM),  we  need 
evaluate  the  integrals 

,  f(,VrV-'tR 

I(,)(r')=f— ^ - dS,  f<Z)(r)  =  i*  — 17  ,  /  =  1,2,3  (1.1) 

•  A 


where 


R  =  \Jr'2  -2r' -r  +  r?  =aJ(x-x')  +(v-v')  +  h1  (1.2) 

with 

r  =  xx  +  yy,  r '  =  x'x  +  y'y  +  hz  (1.3) 

and  the  region  of  integration  T  being  the  triangle  in  figure  1.1.  The  origin  of  coordinates  is  the 
centroid  of  the  triangle,  and  the  latter  lies  on  the  xv-plane.  Boldface  letters  denote  vectors.  The 
same  letters  in  italics  denote  the  magnitudes  of  these  vectors  while  a  caret  over  a  letter  denotes  a 
unit  vector.  The  vector  r;  is  the  position  vector  to  the  /-the  vertex  of  the  triangle,  as  shown  in 
figure  1.1. 


Figure  1.1:  The  integration  triangle 
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We  note  that 


r(0 


(*-')  =  J 


(r-r,)e-" 
R 


■^=J 


(r-r')c  '** 
R 

-ikR 


■dS-(  r,  -  r ')  J  - 


-ikR 


R 


-dS 


dS 


-ikR 


-ikR 


<r'  -p')|  + M dS = -  *  1 ' R"  )ds  -(r'  -  p')J 


-ikR 


R 


ik  ■ 


R 


-dS 


(1.4) 


where  v  is  the  exterior  unit  nonnal  to  the  side  of  a  triangle,  lying  on  the  plane  of  the  triangle. 
The  last  integral  in  this  expression  (the  Scalar  integral)  also  appears  by  itself  in  the  MoM  and, 
thus,  it  is  one  of  the  integrals  that  needs  to  be  evaluated.  The  integral  before  it  can  be  converted 
to  a  line  integral  around  the  triangle’s  boundary  using  a  standard  identity  ([2],  p.  503) 


r(7) 


=  J  (r  ^  dS  =  ^  J  ^  ^  )dS  ~  fc  ~  P')  J 

XV.  I*  *V  j. 


-ikR 


R 


-dS 


(1.5) 


We  can  proceed  one  step  further  and  measure  all  distances  in  wavelengths.  We  note  that 

fds'^Afds,  jdS'^A2jdS  (1.6) 

dT  dT  T'  T 


when  we  move  from  measuring  length  in  meters  to  length  in  wavelengths.  With  this  in  mind,  and 
dividing  all  distances  by  wavelength,  we  get  from  ( 1 .5) 


i(,V) 


-i2nR 


-UnR 

)ds-A2(rl-p')\^-dS 


(1.7) 


where  all  distances,  including  triangle  dimensions,  are  measured  in  wavelengths.  We  note  that 
both  integrals  are  independent  of  the  index  /,  and  that,  although  the  first  one  is  a  vector  integral, 
it  is  in  essence  a  vector  sum  of  three  Scalar  integrals,  each  defined  over  a  side  of  the  triangle. 

The  objective  of  this  report  is  to  calculate  the  integrals  in  (1.7)  to  a  pre-assigned  number  of 
significant  digits.  Neither  of  these  integrals  can  be  evaluated  analytically.  The  integrands, 
however,  can  be  replaced  by  functions  that  we  know  how  to  integrate  and  that  can  approximate 
the  actual  integrands  arbitrarily  closely.  The  obvious  way  to  proceed  is  to  expand  the  exponential 
function  in  a  Maclaurin  series  about  R  =  0  and  truncate  after  the  required  number  of  terms.  This 
works  well  when  R  is  small,  so  that  the  number  of  terms  for  agreement  to  the  required  number  of 
significant  digits  is  also  small.  This  translates  to  the  observation  point  being  near  the  integration 
triangle.  The  observation  point,  however,  can  be  any  point  in  space,  except  for  points  on  the 
three  sides  of  the  triangle.  Thus,  as  the  observation  point  recedes  from  the  triangle,  more  and 
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more  terms  are  required  in  the  Maclaurin  sum  and  there  comes  a  point  when  this  method 
becomes  counter-productive  due  to  the  large  amount  of  time  required  to  compute  this  sum.  We 
attempt  to  get  around  this  difficulty  by  multiplying  and  dividing  ( 1 .7)  by  a  common  factor;  we 
thus  write 


-iljur' 


I</,(r')  =  -A2 


i2n 


dT 


—i27u(R—r'} 


i2n{R— r') 


ds-X 2(r,  -p')e12^  ^ - dS  . 


(1-8) 


We  can  show  that  the  argument  of  the  exponential  in  this  case  is  bounded.  The  centroid  of  a 
triangle  divides  each  median  in  a  2:1  ratio.  The  longer  of  the  two  parts  is  the  distance  from  the 
centroid  to  the  vertex  of  the  median.  Thus,  the  largest  distance  from  the  centroid  to  a  point  on  the 
triangle  is  equal  to  two  thirds  the  longest  of  the  three  medians.  But  a  median  is  shorter  than  the 
longer  of  the  two  adjacent  sides;  thus,  the  longest  distance  from  the  centroid  to  the  triangle’s 
perimeter  is  less  than  two  thirds  the  length  of  the  triangle’s  longest  side.  We  combine  this  with 
the  fact  that  the  side  of  a  triangle  (in  this  case  that  formed  by  the  observation,  integration  and 
centroid  points)  is  always  greater  than  the  difference  of  the  other  two  sides  to  obtain  the  bound 

\R~r'\Zr<hm  (1.9) 


where  /max  is  the  length  of  the  longest  side  of  the  triangle.  Though  this  is  not  a  strict  bound,  it  is 

a  bound  that  holds  for  all  observation  points;  thus,  no  matter  how  far  away  the  observation  point 
is  from  the  integration  triangle,  the  number  of  terms  in  the  expansion  will  be  the  same  as  for  a 
point  near  the  triangle.  Moreover,  for  an  actual  grid,  we  can  search  among  all  triangles  for  the 
longest  side  and  use  that  value  in  (1.9).  In  this  way,  we  do  not  have  to  test  all  triangles 
separately,  conserving  a  lot  of  compute  time. 

We  proceed  now  to  expand  ( 1 .8)  in  a  Maclaurin  series  and  keep  the  first  N  terms 


I^(r')  =  T2e 


-i2nr' 


T 

n= 0 


\(-nx)n  f0 R-r'Y 


n\ 


-IviR-rJds-X'h- p')e^Z^4 

Qf  n= 0  ^  •  x 


R 


</S.(1.10) 


In  determining  the  number  of  terms  N,  we  use  Taylor’s  Theorem  with  a  Remainder  ([3],  p.  1 13). 
If  we  require  a  number  M  of  correct  significant  digits,  then  we  proceed  as  follows  to  determine 
the  number  of  tenns.  The  series  for  the  sine  is 


sin  2k(R  -/*') 
2  K(R-r') 


^(-l)"[2;r(R -/)]"" 

h  (2«  +  l)! 


+  R,(N,2x(R-r')) 


(1.11) 
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where  i?s  is  the  remainder  of  the  series  and  is  bounded  by  the  first  omitted  term 


< 


i  fi  ~\2N 

2  7c\R-r\  \ 

I  Imax  J 


f  An,  ' 2N 


< 


[2N  + 1)!  (2N  +  1)! 


(1.12) 


We  can  then  write  for  the  relative  error 


sin[2;r(f?- 

r')]  ^(-iJ-f^^-rOf 

2  Tr(R-r')  ^  (2n  +  l)\ 

sin[2;r(f?-r')] 

A 

1 

21 

<N 

R,(N,2n(R-r')) 


A, 


IN 


A, 


.IN 


sin[2;r(f?-r')] 

-< - 1 

(27V  +  1)! 

sin \l.n  7?  —  r' 

max  _ 

sin 

(4  7T  i  ' 

A 

1 

21 

_ a 

A 

1 

<N 

[2N  +  1)! 

a  max 

v  7 

ax 

4  n  , 

(1.13) 


The  number  of  terms  is  determined  by  requiring  that  the  relative  error  is  smaller  than  5  divided 
by  the  number  10  raised  to  the  number  of  significant  digits  Mplus  one,  or 


rs  max 

4 - 4—  <  5  -10 

(2JV  +  1)! 


sin 


4  n 


A, 


A 

) 


(1.14) 


Once  the  triangle’s  dimensions  are  known,  we  can  solve  this  expression  for  N. 

In  this  discussion  we  have  assumed  that  the  argument  of  the  sine  function  is  small;  in  practice, 
the  longest  side  of  a  triangle  does  not  exceed  one  tenth  of  a  wavelength  and,  hence,  the  largest 
argument  is  approximately  equal  to  0.4tt/3  radians  or  24  deg.  This  means  that  the  trigonometric 
functions  in  the  original  integrals  hardly  exhibit  an  oscillatory  behavior. 

The  argument  for  the  cosine  function  is  slightly  more  complicated  and  is  driven  by  the  fact 
that  the  integral  of  the  leading  tenn  of  the  first  sum  in  (1.10)  is  equal  to  zero.  This  can  be  seen 
most  clearly  in  (1.4)  where  the  leading  term  in  the  expansion  is  a  constant  (1)  and,  hence,  its 
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surface  gradient  is  zero.  Thus,  the  first  contributor  to  the  integral  is  the  second  term  in  the  cosine 
expansion.  We  then  proceed  to  re-write  the  relevant  integral  in  the  fonn 

I  =  J  i/{c°s[^(i?-r')]-l|(i5  .  (1.15) 

ST 


We  expand  now  the  new  integrand  into  a  Maclaurin  series 


r  (— 1 )"  f 2zr(i?  —  ^')T”  r  -,2  -A  (-1)T2;t(.K -/'')] 

cos[2^(R-/)]-1  =  £V  1  - —  =  -[MR-r,)]IJ  J  (1-16) 


n—\ 


(H 


n= 0 


(2/7  +  2)! 


We  consider  now  the  termination  of  the  last  power  series  after  N  terms 


l-cos[2.(^-,')]  _ g tTlMATir  +Rc(n,2k(R- r')) . 


\_2ir(R-r'j\ 


(2/7  +  2)! 


(1.17) 


We  plot  the  left-hand  side  of  this  equation  in  figure  1.2.  We  see  that,  in  the  range  of  interest,  it  is 
strictly  decreasing  as  its  argument  increases. 
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The  remainder  in  (1 . 17)  is  bounded  by 


R(N.2K(R-r')) 


J 

k\R-r'\ 

1  Imax  , 

rj 

f 

^  ^max 

(2N  +  2)! 

(IN  +  2)! 

(1.18) 


Taking  to  account  the  monotonicity  of  the  function,  we  can  write  for  the  relative  error 

f  A  n 


,2  N 


Rc{N,2K(R-r')) 


-L 


< 


max 

v  j  J 


1-cos  2 n(R-r'Y^ 

i  ( 47r  i  1 

1  — cos  — / 

[2n(R  -r')]2 

(2N  +  2)! 

‘'max 

V  ^  J 

f  An  V 

^max 

V  4  7 

(1.19) 


and  we  set  the  criterion  for  determining  N  to  be 


A, 


v  2  N 


(2A  +  2)! 


<5  10 


-(M+l) 


1-cos 


A, 


V 


^  An  ^ 

3fmax 

y 


(1.20) 


Using  criteria  (1.14)  and  (1.20)  for  the  sine  and  cosine,  and  assuming  that  the  largest  side  of  a 
triangle  is  less  or  equal  to  one  tenth  of  a  wavelength,  we  arrive  at  table  1.1.  This  table  gives  the 
number  of  terms  required  to  guarantee  a  certain  number  of  significant  digits  (SD).  We  see  that, 
for  single  precision  (7  SD),  we  need  a  total  of  8  terms  while,  for  double  (15  SD),  we  need  a  total 
of  14.  We  see  that  only  on  one  occasion  (6  SD)  does  the  number  of  terms  between  sine  and 
cosine  differ.  Finally,  the  criteria  we  use  are  universal,  i.e.,  they  are  independent  of  the  location 
of  the  observation  point  and  the  shape  of  the  triangle.  Thus,  once  we  decide  on  the  number  of 
significant  digits,  we  can  calculate  the  number  of  terms  once  and  for  all  provided  we  know  the 
largest  side  among  all  triangles  in  a  grid. 


9 


NAWCADPAX/TR-2008/227 


Table  1.1:  Number  of  terms  required  to  guarantee  a  given  number  of  significant  digits  (SD). 
Largest  side  of  triangle  less  than  or  equal  to  a  tenth  of  a  wavelength. 


Number  of  SD 

Sine 

Cosine 

4 

3 

3 

5 

3 

3 

6 

4 

3 

7 

4 

4 

8 

4 

4 

9 

5 

5 

10 

5 

5 

11 

5 

5 

12 

6 

6 

13 

6 

6 

14 

6 

6 

15 

7 

7 

As  mentioned  above,  using  the  same  identity  that  allowed  us  to  convert  the  first  integral  to  a 
line  integral,  we  can  show  that  the  zeroth-order  term  in  the  line  integral  of  (1.10)  does  not 
contribute.  Down-shifting  the  index  by  one  we  get 


N-2 


(-n^y 


-.(-&)■  .(R-rJ 


(n  + 1)!  L  n\  J 


dT 


R 


dS 


h  {n  + 1)! 


w(r')-^(r-p  '  =  1.2.3 


(1.21) 


where 


!(.,(>•')  =  j  0(R  -  rT'  *  .  *<„)(■•')  =  (»  +  l)j {R  r)  dS . 

dT  T  K 


(1.22) 


It  is  these  two  integrals  that  we  must  compute  analytically.  The  first  is  a  line  integral  while  the 
second  is  a  surface  integral.  In  Appendix  A  we  employ  the  Gordon-Bilow  transformation  to 
convert  the  surface  integral  to  a  line  integral  [4],  [5].  From  (A.l)  and  (A.42),  we  have  that 


sj 


P~ 


(M-A 


-ds 


(1.23) 
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where  s  j  is  the  y'-th  side  of  the  triangle  and  vj  is  the  exterior  unit  normal  to  it,  lying  on  the 

triangle’s  plane,  while  p  is  the  distance  between  the  integration  point  and  the  projection  of  the 
observation  point  on  the  triangle’s  plane.  In  this  notation 

lw(r-)  =  ±vJUJpr^-rf'ds.  (1.24) 

We  evaluate  these  integrals  in  the  next  section. 

2.  EVALUATION  OF  LINE  INTEGRALS 
We  begin  by  writing  (1.23)  in  the  form 

vr'>=i|v(r/.,-p'KJ(o  (2.i) 


where 


k,A)=J 


(7p2+A2  -r) 


-ds ,  n  =  0,1,2,.... 


P 


(2.2) 


From  this 


Ku  =  1 


\[p2  +/r  - 1*| 


ds 


P~ 


K\j  =  Sj  -  2r'K0  j 

Kn+lJ=Vn_lJ-2r'KnJ-p'2Kn_lJ,  n  =  1,2,. 


where 


F, 


,j(r')  =  \Up2  +h2  -r'\  ds  ,  n  =  0,1,2,. 


(2.3) 

(2.4) 

(2.5) 


(2.6) 


Comparing  the  last  integral  to  that  in  (1.24),  we  find  that  we  can  write 

t.,  (r')  =  I>/.j(r')-  (2-7) 

y=i 
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Before  proceeding  with  the  evaluation  of  the  integrals  in  (2.6),  we  summarize  some  essential 
definitions 

r  =  r/+l  +  stj ,  j  =  1,2,3  (2.8) 

with  the  indices  running  cyclically.  We  also  write 

r'  =  x'x  +  y'y  +  hz  =  p'  +  hz  ,  r,  =  x,x  +  yty .  (2.9) 


and 


ay+1  =  rJ+1  -  p' ,  T  =  s  +  tj  •  a  ,+1  ,  A.  =  t}  x  (it  x  a  .+1 )  ,  p  =  |r  -  p'|  =  |a  .+I  +  stj \ .  (2.10) 


From  these  we  can  write 


P2  =s2+  2 ij  ■  a J+ls  +  |a  ,+1 1~  =  (s  +  ij  ■  a /+1 )"  +  |a ,+1 1  -  (fy  •  a  .+1 )  =  (j  +  •  a  .+1 )  +  =  r2  +  A  2 

(2.11) 


We  also  define 

by  =r/  -r';  Bj=Jaj2  +  h2  ,  j  =  1,2,3  (2.12) 

The  geometric  significance  of  these  quantities  is  shown  in  figures.  2.1  and  2.2. 
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Vertex  j 


Figure  2.1 :  Geometrical  meaning  of  the  quantities  in  Equation  (2.10) 

2.1  EVALUATION  OF  Vn  ] 


For  ij  ■  b  /+2  >  ij  ■  b  /:1  >  0 ,  we  make  the  transformation 


t  =  Bj  tan# ,  =  #,  sec"  ### 


and  write  for  (2.6) 


M  = 


f  (Jr2+5.2 -/V 

t+1 

dz  -  B 

K  J  J 

i  Msec#-/) 

j 

J  V  ]  J 

/ ~  .  \ 

n  ti 

-I 


n= 0 


(n  +  ^  m+ 1  /_  A»+1-m 

m !(/?  +  !-  m ) !  ^  7 


J  secm+2  6d0 . 


V/+ 1 


(2.13) 


(2.14) 
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Vertex  j 


Let 


Figure  2.2:  Local  rectangular  coordinates  and  various  vectors 


v  =B 

m  ,J  J 


n+\ 


J  sec  m+2dd9,  m  =  0,1,2,... 


Wi 


(2.15) 


Then 


vo  j=Bj 


sec  6d6  =  B  tan# 


while 


arctan 

r  t  ■  b  ^ 

lj  uj+ 2 

B 

V  J  J 

arctan 

B 

V  J  J 

=  t..bJ+2-tj-bj+l=sj 


(2.16) 


arctan  '’"b'+2 


=  BJ 


J  sec3  6d0  =  B2  sec  6  tan  6 


arctan  'j'bj+l 


+  A:‘ 


f 

7 

arctan 

fj  -b,+2 

Ih2  +  A2 

V  J 

7 

f 

arctan 

h  •  V i 

,  , 

Jh2  +  A2 

v7 

V  2 

y 
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=  B 


1  + 


arctanj 

,  J 

r 

arctan 

J 

{ 7A+i) 

l  *j  ) 

(h  • 

\2 

•vh 

B 


1/2 


V  BJ  J 


1  + 


ft-M2 

B2 


1/2 


v  bj  , 


—  v, 


=  b./+2  |  (O  '  b  y+2  )  -  |b7+l  |  (h  '  b/+l  )  ~  Vl,7  +  B)  J  SeC  0d0 


7'V  l 

B, 


1 J 


+A  j 


so  that 


vu  = 


+  5; 


=+!K.|vb+2)-lbMlft'b+.) 


In 


1  + 


h  -b./+2 

V  BJ  J 


+ 


h  -b/+2 

V  B>  J 


-In,  1  + 


ft  b  V  f 

L]  uj+ 1 


v  »j  j 


+ 


(j  'b/+i 


V  BJ  J 


=  -< 

2 


b7+2|ft-b7+2)-|b7+l|ft-b7+l)  +  57 


In 


b7+2  +*j  '  b  /+2 

B 


-In 


by+i  +tJ  ~  b/+i 

B 


by •  2  ft-b/-0-b/  i  ('Vb/-i)  +  /?/2 1,1 


b./+2 

+  tj 

•b./+2 

by+i 

+  tj 

‘V 

sec  6d6 


(2.17) 


(2.18) 
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In  general, 


v  =B 

J 


m+ 1 


J  sec”'+2  6d6  -  B'n+l  sec'"  0tan0 


arctan 

T-Om) 

B 

arctan 

(+■>*,) 

B 

v  2  y 

-/ nBjm+1  |  sec'"  0  tan2  0d0  =  /?/""  sec'"  0  tan  0 


arctan 


arctan 


f  t  ■  b  ^ 

'j  uj+ 2 


„  J  m 

h  •*>/+■ 

V  *J  J 


(v  —  B  2v  ,  . ) 

\  miJ  J  m—2 ,j  J 

(2.19) 


or 


B 


m+l 


"uj  m  + 1 


1  + 


^•b,+2V 

v  ^  y 


-\m/ 2 


^Vb/+ 

v  y 


1 

m  + 1 


by+2  (*/ 


1  + 


m 


v  y 


-\m/2 


'tj  -b./n  ^ 
v  y 


+ - -B,2vm_2i,  m  =  2,3,... . 

m  +  l 


)l  \m  / * 

~  |by+i  ('/  + 

This  completes  the  evaluation  of  (2.15).  From  (2.14)  and  (2.15) 

A2+1 

f«j=Zi  ././  x. k-^ ; 


m= 0 


(  77  +  1 ) 

!  1 

m!( 

77  +  1- 

-777 

U 

V„!y.,  n  =  0,1,2, 


+  - 


mB2 
m  + 1 


K  ~  • 

m— 2,7 


(2.20) 


(2.21) 


We  examine  now  the  case  where  £.  •  b/+l  <  £.  ■  b;+2  <  0  .  From  (2.6) 


^  = 


*j  'b;+i 


=  |  +  ZT2  -r'j  dr  =  J  +  ~bJ  —r'  I  dz 


n+ 1 


/;-b 


ri-ri 

-I 


(7?  + 1)! 


/'uy+2| 
arctan  I 


777 !  ( n  +  1  -  in ) ! 


j)  m+ 1  /  Aw+1  m  f  m+2  /n  i 

Bj  [-r  j  I  sec  OdO . 


(2.22) 
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Proceeding  as  above, 


i+ 1 

B, 


D  n 

V  —  B 

™,J  J 


J  sec  m+20d6,  m  =  0,1,2,... 


'i'uj+2 | 


(2.23) 


'  0'by+i  | 

arctan 

f  tj '  b/+. 

\  ) 

B 

voj  = 


fF  J  sec2  Odd  =  B/  tan# 


arctan 


tJ  -V  2 

5 


<=  VVi  'b/+2 


(2.24) 


vu  = 


b/+1  ■  b  /+i  -  b  /+2  tj  ■  bj+2  +Bj  In 


bi+1  +  tj  -by+1 
b/+2  +  tj  ■  b  /+2 


(2.25) 


v  = - 

m,/  ^ 

in  + 1 


V  h  'by+i  -  by+2  h  '  b;+2 


^  B  2  o  -j 

+ - B,vm_2,,  m  =  2,3,... 


m  + 1 


(2.26) 


For  the  case, 


tj  •  b,+i  <  0  ,  tj  ■  bj+2  >  0 


(2.27) 


we  write  from  (2.6) 


v*.j  =  j  (Jt2+Bj2  ~r')  dT  =  J  (^^T^+Bj2-r,)  dT+  I  (ylr2+Bj2  ~r') 


t  +  Br  -r'\  dr 


-  1 


t  4-  B,  -r'\  dr  + 


n+ 1 


i  (4 


t  +  B .  -/ 


K+l 


W+l 

I 

m=0 


(n  + 1)! 


m !  ( n  +  !-/«)! 


B  m+l  \ 


J  sec"i+2  OdQ+  J  sec m+1 6dQ 


.  (2.28) 
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In  this  case,  (2.15)  becomes 


D  /W+l  I 

v  =  B  < 

m,j  j  ’ 


l  Vb/4 

1  \ 

•i| 

a  rc  tan 

fOA+a) 

I  ^ 

/ 

l  Bi  ) 

v  j  /  \  j  y 

J  sec m+20d0+  |  sec"1+2  OdO 


m  =  0, 1, 2, . 


(2.29) 


so  that 


vo  ,j=BJ< 


J  sec n6d0+  J  sec  Odd 


r=  tj  ■  b  /+l  +tj  •  b  /+2  =  sy 


(2.30) 


\J  2  < 


by+i  1 1 tj  ■  by+ 1 1  +  |by+2  |  (o  •  by+2  )  +  ln 


(|b2'+l  |  +  jt/'  ‘  bt+l  |)(|bi+2  \+tj  '  by+2  ) 


B 


(2.31) 


v  = - 

m.  1  < 

m  + 1 


7+1 


h  *  by+l 


+  b./+2 


tVb/+2) 


+  - 


111 

m  + 1 


Bj\- 2J  , 


m  =  2,3,... 


(2.32) 


This  concludes  the  evaluation  of  (2.6).  The  iteration  formulas  (2.20),  (2.26)  and  (2.32)  do  not 
converge  everywhere.  In  Appendix  B,  we  show  that,  for  the  iteration  to  be  stable,  we  need 

0  <  Bj .  <  1 . 


We  also  mention  that  there  are  other  ways  of  evaluating  these  integrals  but  they  do  not  appear 
to  have  an  advantage  over  the  present  one. 

2.2.  EVALUATION  OF  K0J 


We  turn  to  the  evaluation  of  (2.3)  and  write 

(/r  +h2y  —  \h 


Ko  j  = 


J' 


P~ 


ds  =  \ 


ds 

-  1 

p  ds 

(p2+A2) 

1  J 
p+H  1 

U+w 

(2.33) 
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We  make  the  transformation 

R  =(s  +  tj  ■  a/+1  )  +  Bj  => 5  +  fy  •  a j+l 


ds  = 


±RdR 


(2.34) 


and  consider  observation  points  in  the  region  f .  •  b  /+2  >  f .  •  b  /+l  >  0 .  In  this  case,  the  appropriate 
sign  in  (2.34)  is  the  positive  one.  From  the  last  two  statements 


K 


0,7 


j*|V2|  RdR 
■'M  (R  +  \h\)yjR2  -  Bj2 


=  A 


-h 


(2.35) 


We  let 


R  =  Bf  cosh(w) ,  Ji?  =  Ry  sinh(n)Jw 
and  write 


ft 


=  In 


=  In 


b/+2 

dR 

cosh 

-J 

V  i 

\ 

Ir2- 

-B2 

j 

cosh-1 

b7+2 

W 

b7+2 

2  9 

~Bj 

B 

J 

b7+2 

+\l 

b7+2 

2  9 

-Bj 

7+1 1 

B: 


du  =  cosh  1 


V  *,  ) 


-cosh  1 


v  *,  j 


\  f 

-In 


b7+l  + 


'7+1 


-5 


5, 


'7+1 


b,,  +JbJ  -5, 


7+1  7 


=  In 


b7+2 

+  ij '  b7+2 

b7+l 

+  ij '  b7+i 

For  the  second  integral  in  (2.35),  we  write 


r-Jt  ^ 

II 

"  — 

(R-h\)dR 

J  1  'ii,i  (*+H)A2-A 

(*2-H=)A:-A 

-Ih'Y'  RdR 

b7'+2 

-i*r  i 

M 

dR  , 

(R2-h2)^R2-B/  21  " 

(2.36) 


(2.37) 


(2.38) 


19 


NAWCADPAX/TR-2008/227 


In  order  to  evaluate  I2l ,  we  let 
u2  =  R2-B 2 


(2.39) 


du 

u 1  +  A 2 


(2.40) 


We  now  make  the  transformation 


u  =  A.  tan i// 


(2.41) 


to  get 


21 


I  y  J-* 

=1*1  J 


K;  ‘b 


,/'Dy'+i 


du  _  h 
u2  +  A2  A . 


J  J 

J  dy/  = 


4 


tan 


■  b  ^ 

v  uy+2 

A,. 


-  tan 


*4  b 

Lj  ui+ 1 

A,. 


(2.42) 


For  the  integral  /22  in  (2.38),  we  let 


/? 


v  = 


(R2-R2)  =  R2  =>  /?  =  £ . 


•n/v2  -1 


(2.43) 


so  that 


4^1- 


dR  = 


Vv'’  -1 


v2-l 


dv  =  — 


5. 


/  ,  \  3/2 

(v  -  0 


dv . 


(2.44) 
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We  can  then  write  for  I21  in  (2.38) 


4=W2  J 


dR 


bj+ 2 
IVb7+2 


l\(R2-\h\2U 


R2  -Bf 


=  ~\ht  J 


dv 


m-2 


\tj-b 


J-°J+ 1 


(2.45) 


We  let  next 


A. 


W  =  -r—V 


\k\ 

and  write 


(2.46) 


\tr  b 


j’°j+ 2 


4=-i*r  i 


h  dw 


(4v)~  + 4 


A 


w2  +1 


\h\ 

f 

1  1 

tan-1 

A. 

j 

V 

4  Vi 


-  tan- 


f  a  k  ^ 

4  °j+2 


Hh  -V 

From  (2.38),  (2.42)  and  (2.47),  we  have 


VI h  h  'bj+2\j 


(2.47) 


I 2  1 21  ‘22 


— 

r  ~  > 

— 

\h\ 

tan-1 

t  •  b 

lj  4+ 2 

-tan  1 

t  -b 

V  4+1 

4 

4 

A. 

J 

l  '  7 

K.  J  J 

\h\ 

+  — 

tan-1 

A. 

j 

V 

4  4+2 
I/2I  b, 


-  tan 


f  a  b  ' 
4  4+i 

vi^l|4Vi7 


(2.48) 


From  this,  (2.35)  and  (2.37),  we  can  write  that 


Koj=I\-I2=  In 


4+2 

+ 

h  ■  4+2 

4+> 

+ 

h  -4+i 

A, 


tan 


b 

Lj  4+2 

A. 


-  tan 


ni  b  ^ 
Lj  4+ 1 

A, 


lAl 

1  1 

tan-1 

7 

V 

4  4+ 2 

\h\  t,.  b  ;j+ 


-  tan 


vi/2I|4Vi; 


tj  ‘  4+2  >  t]  ■  by+i  >  ° . 


(2.49) 
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As  we  pointed  out  above,  this  formula  is  valid  when  the  integration  in  (2.33)  is  over  an 
interval  of  positive  numbers.  This  is  made  explicit  through  the  inequality  in  (2.49).  Basically, 
and  with  respect  to  the y'-th  side  of  the  triangle,  we  divide  all  of  space  into  three  sectors.  They  are 
formed  by  attaching  an  infinite  plane  at  each  end  of  the  y-th  side  and  normal  to  it.  The  inequality 
in  (2.49)  holds  in  one  of  the  two  sectors  that  do  not  contain  the  y-th  side.  In  the  other  sector,  the 
directions  of  the  inequality  are  reversed.  In  this  case,  and  by  (2.34),  we  must  introduce  a  negative 
sign  in  (2.35)  or,  simply,  reverse  the  direction  of  integration;  thus,  from  (2.49)  we  have 


K0j  =  In 


V 

+ 

ij  ‘V 

by+2 

+ 

ij '  b  /+2 

M 

A,. 


tan 


f  A  b  ^ 

Aj  °j+2 


AJ 


tan  1 


ni  -b  IA 

li  "j+ 2 

A,. 


-tan  1 


ni  -b  IA 

vj+ i 

A, 


vri|0  'b/+2|y 


-tan  1 


f  A  b  ^ 
Aj  Dy+i 


\h\tj-b 


J«\J 


tj  •  b./+i  <  tj  •  b/+2  <  0  • 


(2.50) 


The  third  case  is  when  the  observation  point  is  in  the  space  between  the  two  planes.  The 
corresponding  inequalities  are 


tj  ■  b/+,  <  0  ,  tj  -bi+2>0 


(2.51) 


which  necessitate  splitting  the  interval  of  integration  in  two.  The  resulting  expression  is 
RdR  r b  j + 2  RdR 


r  b. 


d+nU 


:  + 


r2-b 


2  J  B 


l 

JB: 


R--B 


=  In 


by+1  +  tj  •  b  /+1 

B 


+  ln 


nb  +  f  b  ^ 

uj+ 2  ^  Lj  uj+ 2 

B; 


---tan1 

Aj 

\h\  j 
— ■ — -  tan 
A, 


h  -b/+. 

As 


h  '  by+2 
A: 


M 

aj 

W 

aj 


tan  1 


tan 


f  A  b  ^ 
Aj  P,/+l 

VH|0  ‘  bt+i  |  J 

f  A  |b  I  ^ 

y\h\  tj  ■  by+2  |  j 


n 

~2 

n 

~2 


=  In 


(K.i|+l'Vv,|)(KAI'Vb^l) 


M 

A, 


cot 


B 


f  A  b  ^ 
Aj  Dy+i 


M 

A 


tan 


ni  b  |A 

2  vj+ i 
A.. 


+  tan 


At  -b  ^ 

2  Uy+ 2 

A, 


\h\ij-b 


+  cot- 


J*u 


f  A  b  ^ 

Ai  Pt+2 

yM  ij  ' by+2 1 y 


=  In 


(|*Vi  |  +  I  ij  •  by+i  |)  (|b7+2 1  +  \ij  •  by+2 1) 


B 


A 


tan 


ni  -b  IA 

A, 


+  tan 


ni  -b  IA 
2  2+2 

A, 
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A, 


tan 


f i,  ;  ,  A 

h  tj  •  b /+1 

V  a>K\  J 


+  tan  1 


j  Q  '4  + 2 

v  4K-I  ; 


tj  •  b  /+i  <  °  ,  tj  •  b/+2  >  0 , 


We  revisit  (2.49)  and  combine  the  first  and  third  arctangents  according  to  the  fonnula 

f  „  ,  ,A 

tan  1 

tan-1  (x)  +  tan-1  ( v)  =  \ 


x  +  y 


n  +  tan 


v1-^; 

f  -»/•  i  t .  ^\ 

-1 


*  +  y 

l-xyj 


xv  <1 

x  >  0  ,  xv  >  1 


Since 


1 J  * by+2  4  4-2 


'i+2 


AJ  * b /-2 


>1 


we  use  the  lower  fonnula  in  (2.53)  and  write 


tan 


nt  b 
v  uy+z 

A 


+  tan  1 


f  ^  b  ^ 
4  P./+2 

Vl/2|  tj  ■ bj+2\J 


=  7i  +  tan  1 


l,  '  4+2  ^  4 

b  y'+2 

4 

A 

h  *  b2+2 

h 

1-L 

'y+2 

A 

=  n  +  tan 


h\ - +  — - 

Aj  tj  •  b  /+2 


h-  b 


A  2 


=  n  -  tan 


I,  I  ^ j  '  b7+2  4+2 

« - b 


4 


'4+ 2 


2  2 
b./+2  “|* 


(4|  +  |b./+2|) 


=  n  -  tan- 


=  #  -  tan- 


=  n  -  tan- 


A 

h  b/2 

2  , 

+4 

by+2  (i 

A  \1  IL 

i2  U 

HI 

Jay'+2 

(W) 

i++/l 

4+2 

1 

2  i* 

L  I2  l 

‘y  t+2  y+2 


KH  +  M  ^ 


4  4  '  4+2 1  14  •  by+2 1 


+  |by+2 1) 


(4  +  by+2  ) 


=  it  -  tan  1 


A  1 

'  h  + 

4+2 

)+42 

x4. 

J 

h  -4+2 

(2.52) 


(2.53) 


(2.54) 


(2.55) 
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Similarly, 


tan 


v  Aj  J 


+  tan- 


(  A  b  ^ 

4  °J+ 1 


=  n  —  tan 


h  | 

4  + 

V 

)+f 

Aj 

h  ~b.M 

Substitution  of  the  last  two  expressions  in  (2.49)  yields 


K0J  =  In 


by+2 

+ 

ij  -b  j+2 

p 

\h\ 

J _ L 

V 

+ 

h  -V 

j 

1 

4 

h  | 

4  + 

V 

by+2 

)+f 

4 

h  ’  b/+2 

-tan 


h(h  + 

Vi 

)+a; 

4 

h 

tan  1 


4  h  '  b  /+2  >  I/ '  b  /+l  >  0 . 


As  Aj  — >  0  ,  we  use  the  well  known  expansion 


tan  1  (x) 


2 


*  (— i)A 

S(2k  +  l)x2*+1  ’ 


to  obtain  an  estimate  for  A.  so  that  only  the  first  term  in  the  series 
specified  number  of  significant  digits.  We  thus  set 


4 

h '  b/+2 

h  | 

4  + 
\ 

by+2 

)+f 

V3-1(TM 


and  by  means  of 


4 

ij '  b  ./+2 

4 

4 

by+2 

h  1 

4  + 

by+2 

)+A 

A| 

4  + 

by+2  ) 

(2.56) 

(2.57) 

(2.58) 

will  be  required  for  a 

(2.59) 

(2.60) 
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we  obtain  the  estimate 


h  1 

4  + 

k 

b^) 

h 

b 

j+2 

A,  <  V3  •  10 


For  these  values  of  A. ,  we  replace  (2.57)  by 


b,+2 

+ 

h  -b/+2 

p 

w 

n 

4  h  V 1 

V 

+ 

2 

2  A) 

^  +  Vi  )+42 

n 

- + 


4 

'  b  /+2 

A| 

'  h  + 

4+2 

=  In 


4+2 

+ 

tj  -b/+2 

p 

-  A  < 

'4+2 

4+  i 

+ 

h  '4+i 

2 

A  | 

l  A  +  4+2  ) 

l+42 

'4+i 

V 

A 

A  + 

4+i 

)+4 

0'by+2  >  tj  ■  by+i  >  0 . 


Moreover,  if  A  =  0 ,  then 


4>,y  =  ln 


(2.62) 


4+2 

+ 

tj  -b,+2 

V 

+ 

h  '  4+i 

,  A  =  0,  tj  ■  bj+2  >  tj  ■  by+1  >  0 . 
Similarly,  in  place  of  (C.  18)  we  write 


-tan- 


by +2 

+ 

‘  4+2 

p 

A 

4+l 

+ 

h  4+1 

J 

A  1 

A  + 

4+2 

)+A 

4 

4b/+2 

*1 

(a  + 

4+i 

1 

(N 

+ 

4 

*j  '  4+1 

tan 


,  i.  •  b  /+l  <  £,  •  b  /+2  <  0 


and,  if 


A  1 

r  A  + 

k 

M 

h 

b 

/+ 1 

(2.61) 


(2.63) 


(2.64) 


(2.65) 
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then 


*<u=-h 


b;+2 

+ 

4 

*  by+2 

7 

+  \h\- 

h '  b/+2 

h  'hM 

V 

+ 

0 

Z+i 

, 

h 

h  +  by+2  ( 

+  42  h 

+42 

(2.66) 


while,  if  h  =  0 , 


by+2 

+ 

ij  b/+2 

V 

+ 

ZV 

,  A  =  0,  £,.  •  b  /+l  <  ij  •  b  /+2  <  0 . 


(2.67) 


In  (2.52),  we  recall  that  if  h  =  0,  then  A  .  >  0 .  Also,  if  A  .  =  0 ,  then  \h\  >  0.  Otherwise,  the 

observation  point  would  lie  on  the  j- th  side  of  the  triangle.  In  this  expression,  we  combine 
arctangents  according  to  the  formula 


tan  1  (x)  -  tan  1  ( v)  =  tan 


x-y 

1  +  xy 


XV  >  -1  . 


(2.68) 


In  place  of  (2.52)  we  write 


K0j=^ 


(|by-+1|  +  |^  *  by+1 1 )  ( |b y+2  |  +  | ^7  *  ^ y+2  | ) 


h 

+ - 

4 

=  In 

\h\ 

~4 

=  In 


tan 


h  Vb,+i 

V  Ai  by+i  ) 


B 


+  tan 


A: 


tan  1 


h  -  Vi 
A, 


\  f 

+  tan-1 


) 


h  h2 

A. 


h  tj  •  b /+2 
V  AjK* 


(|b/+i|  +  |4by+i|)(|bM]  +  |4b.+2|) 


B 


h 

z 


tan- 


tj *  by+i 

A, 


\  ( 

i 


-tan 


h  tj- by+1 

V  4'K+1  J 


tan 


trhn-2 

A, 


\  ( 
-l 


-tan 


h  tj- bv+2 

V  4  |by+2  ; 


(K+i|  +  |Zby+i|)(M  +  |Zby+2|) 


B 


-  —  tan 

A, 


tj '  b; ,]  htj-  b.+1 


4  4 

by+i 

i+ 

0  -by+i 

|/r| 

0  ‘by+i 

4  4  bt+i 
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--tan’1 

4 


tj  '  b  /  -2  k  44+2 


4  4 

4+2 

1+ 

44+2 

41 

^  -4+ 2 

4  4  4+2 


=  In 


(|V|+|Vbtn|)(|by+2|+|4by+2|) 


B 


-  —  tan 
A, 


-\ 


•4+||4+i|-4l) 


A, 


4 2  4+i  +  444+1 


|2 


h 

- tan 

A, 


|V4+2|(|4+2|-4l) 


A,. 


4 2  4+2  +4  tj  '4+2 


|2  “J 


=  In 


(|4+i|+|44+i|)(M'f|44+2|) 


B 


h 

- tan 

A 

j 


•4+2|(|4+2|-4) 


-  —  tan 
A, 


A, 


•4+i|(|bM]-4) 


42  4+i  +4  4+4-42-; 


42  4+2  +4  4+2 2 -42-*2 


=  In 


(|4+i|+|v4+i|)(|4+2|+|44+2|) 


B 


-  —  tan 

A, 


-1 


44+1  4 


-—tan 
A, 


-1 


tj  -4+2  4 


Aj 2  +41(4+2  +41) 


4" +4l(|4+i|+4 

,  44+, <0,  tj -b;.+2 > 0 . 


We  consider  now  the  case  A  — >  0  .  From 


*  (-1)*  x2i+1  , 

tan  (x)  =  V - ,  jc  <  1 

v  ’  h  2^+1  1 


we  have  that 


tan 


1  (lO'8)  =  10 


1016  1024 
1 - + - 


(2.69) 


(2.70) 


(2.71) 
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Thus,  for  arguments  less  than  or  equal  to  10  ,  the  second  and  subsequent  terms  do  not  contribute 
in  a  machine  that  uses  double  precision.  We  can  then  replace  the  arctangent  by  its  argument.  If 


_ h  ■|V<  4 

A}2  +  \h\(bj+l  +|A|j 


(2.72) 


we  replace  (2.69)  by 


K0J=ln 


>y+i  +  *j 


•Vi|)(M^!‘bAa|) 


tj-bj+l  h 


4  + 


4[|v|+w) 


h  ,  t.  •  b  /+2  A. 

-—tan  '  —A  >  tj'K h-bJ+ 2^° 

4  4  +  4(4+ 2  +  4) 


(2.73) 


while,  if 


tj  -bi+2  Aj 


A:  + 


a|(K+2|+K 


(2.74) 


then  we  replace  (2.69)  by 


K0J=ln 


V 1  +  tj 


•by+1|)(K+2|  +  |Vb.+2[)l  \h 


-  L  tan 

A .. 


tj  ‘by+i  Aj 


A/  + 


_ tj-bj+2  h _ 

Ar  + 1/^|  ( by.+2  +\h\. 


tj-bj+ 1<°,  t,  •  b ,.+2  > 0 . 


j  "j+ 2 


(2.75) 


If  both  (2.72)  and  (2.74)  obtain,  then  we  replace  (2.69)  by 


K0J  =  in 


*7+1  +  t j 


' b /+i |)(|by.+2 1  +  |f,.  -by.+2|) 


tj  ■  b/+!  h 


4  + 


_ tj  'by+2  h _ 

4"  +  ^|(by+2  +|^|) 


4b/+:  -  °,  4b,+2  >0. 


j  "j+ 2 


(2.76) 


28 


NAWCADPAX/TR-2008/227 


If  h  =  0 ,  then  A  j  ^  0  and 


h  =  0,  tj  ■  by+1  <  0  ,  tj  ■  bJ+2  >  0 .  (2.77) 


We  note  that  (2.74)  can  be  used  in  place  of  (2.61),  and  (2.72)  in  place  of  (2.59).  In  connection 
with  this,  we  show  that  these  inequalities  hold  for  both  corners  of  an  edge.  We  first  recall  that 


a  >  b  ,  c  >  0 


a  a  +  c 

—  > - . 

b  b  +  c 


(2.78) 


This  follows  from 


a  >  b  =>  ac>  be  =>  ac  +  ab  >  be  +  ab 
For  the  case  in  (2.57),  consider 

( 


h  -V  2 

V  K'  ‘  b7+i  I J 


b,.+22-(|/*|2+42)  b  12 


b7+1f  -(w2+42 


> 


y+2 


Vi 


a(h  +  c)  >  h(o  +  c) 


from  which 


h 

■b/+2 

h 

'  b/+i 

> 


j+2 


> 


2+1 


h  + 

^  1 

ij  '  b/+2 

i  \ ! 

h  + 

b./+2 

h  + 

b,, 

h  ' b /+i 

h  + 

b./+i 

or 


t]  '  b /+1  ^  tj-kj+2 
h\  +  b  J+l  \h\  +  bj+2 


(2.79) 


(2.80) 


(2.81) 


(2.82) 


If  we  examine  (2.60),  we  see  that,  in  light  of  (2.82),  (2.61)  holds  for  both  vertices  j+l  and  j+2. 
The  same  can  be  said  about  (2.69). 

In  Appendix  C  we  present  an  alternate  way  of  computing  K0  . .  The  rest  of  the  Kn  .  can  be 

computed  from  the  iteration  formulas  (2.4)  and  (2.5).  In  Appendix  C,  we  show  that  the  iteration 
formula  does  not  converge  for  all  positions  of  the  observation  point.  We  discuss  this  issue  in  the 
next  section. 
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3.  COMMON  REGION  OF  VALIDITY  OF  INTEGRAL  EVALUATION 

In  Appendix  B,  we  show  that  the  iteration  formula  (2.5)  converges  provided  (B.14),  namely 


1  —  n'2 

\h\  <  — ,  if  0  <  p'  <  1  and  \h\  >  0 

11  2  11  (3.1) 

0<p'<l  ,  if  \h\  =  0 

holds;  similarly,  (2.20),  (2.26)  and  (2.  32)  converge  provided  (B.19)  is  satisfied: 

0  <  £  <  1 .  (3.2) 

We  also  noted  that  these  two  statements  do  not  imply  one  another.  In  this  section,  we  prove  that 
both  are  satisfied  within  a  sphere  centered  at  the  centroid  of  the  triangle  and  radius  equal  to  0.5. 
The  proof  is  as  follows. 

We  note  that 

r'  =  \] p'2  +h2  <0.5  =>  p'<  0.5,  \h\  <  0.5  (3.3) 

from  which  we  get  that 

|A|< V0.25-P'2  =tVl-V2  <yi-4p’>  +(2p’1f  (3.4) 

so  that  (3.1)  is  satisfied. 

We  also  note  that 

K..|  =  a/V+MZ7  =>  y  =  1,2,3.  (3.5) 

Also,  by  construction,  B .  is  non-negative.  But 

|bi+1|  =  |r'-r.+1|<l  (3.6) 

since  the  observation  point  r' and  any  of  the  three  vertices  r/M  of  the  triangle  are  inside  the 
sphere.  From  (3.5)  and  (3.6)  we  conclude  that  (3.2)  is  satisfied. 

Thus,  our  method  is  guaranteed  to  work  when  the  observation  point  is  less  than  half  a 
wavelength  away  from  the  triangle’s  centroid.  A  question  we  may  ask  in  connection  with  this 
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conclusion  is:  how  does  it  compare  in  speed  with  the  original  Maclaurin  expansion  where  we 
expand  exp (ikR)?  In  table  3.1,  we  present  the  same  information  as  in  table  1.1  except  that  this 
time  we  use  the  unaltered  Maclaurin  series.  The  observation  point  on  the  sphere,  half  a 
wavelength  away  from  the  centroid  of  the  triangle.  If  we  let  it  he  on  the  triangle’s  plane,  then,  in 
the  worst  of  situations,  the  point  on  the  triangle  farthest  from  it  is  the  vertex  with  the  property 
that  (a)  it  corresponds  to  the  smallest  angle  of  the  triangle  and,  (b),  the  straight-line  segment 
connecting  it  to  the  observation  point  passes  through  the  centroid  of  the  triangle.  With  this  in 
mind,  a  bound  on  R  is 


R  <  0.5  +—/ 

3  11 


-I  _2_-lZ 

“  2  +  30~30 


(3.7) 


or 


2  nR<n—.  (3.8) 

15 

We  use  this  as  an  argument  in  (1.14)  and  (1.20)  to  fill  table  3.1  We  see  that,  for  the  observation 
point  on  the  half-wavelength  sphere,  we  require  at  least  twice  as  many  terms  as  in  table  1.1. 
Certainly,  as  the  radius  of  the  sphere  becomes  smaller,  so  is  the  number  of  terms.  From  (3.7)  we 
see  that  R  is  minimized  for  the  sphere  of  zero  radius;  thus,  the  unaltered  method  requires  as  few 
terms  as  the  altered  one  only  when  the  observation  point  is  the  triangle’s  centroid. 

Table  3.1:  Number  of  terms  required  to  guarantee  a  given  number  of  significant  digits  (SD) 
for  the  unaltered  integral.  Largest  side  of  triangle  less  than  or  equal  to  a  tenth  of  a 
wavelength  and  observation  point  is  half  a  wavelength  away  from  centroid. 


Number  of  SD 

Sine 

Cosine 

4 

8 

8 

5 

9 

8 

6 

10 

9 

7 

10 

9 

8 

11 

10 

9 

11 

11 

10 

12 

11 

11 

13 

12 

12 

13 

12 

13 

14 

13 

14 

14 

13 

15 

15 

14 
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4.  VALIDATION 

We  compare  here  results  of  our  method  with  those  of  Khayat  and  Wilton  [6] .  They  compute 
the  second  (Scalar)  integral  in  (1.7),  i.e.,  the  integral  of  the  free-space  Green’s  function 

-ikR 

I(r')  =  \—ds.  (4.1) 

The  integration  triangle  is  shown  in  figure  4.1.  It  is  an  isosceles  right  triangle  whose  equal  legs 
have  length  1  m.  The  wavelength  is  equal  to  10  m.  The  result  of  Rossi  and  Cullen  is  from  [7]  and 
appears  in  table  I  of  [6],  where  it  is  used  as  a  reference  (benchmark).  The  first  Khayat  and 
Wilton  row  refers  to  the  36-point  result  in  table  I  of  [6],  while  the  last  row  comes  from  table  II  in 
[6],  with  a  reported  accuracy  of  14  significant  digits  (SD). 

3:  (-1/3,  2/3,  0) 


Figure  4. 1 :  The  Khayat-Wilton  triangle 

The  results  of  the  comparison  are  shown  in  table  4. 1 .  The  four  observation  points  lie  on  the 
right  angle’s  bisector.  Their  coordinates  are  given  with  respect  to  an  origin  located  at  the  vertex 
of  the  right  angle,  as  in  [6].  According  to  (1.14)  and  (1.20),  the  number  of  terms  required  for 
seven  SD  is  8  (total),  while,  for  15  SD,  it  is  14n  (total). 

In  table  4.2,  we  show  the  difference  between  our  result  for  15  SD  and  that  of  Khayat  and 
Wilton  for  14  SD.  We  have  rounded  our  results  to  14  SD.  There  is  almost  perfect  agreement. 

In  table  4.3,  we  display  the  remaining  points  computed  in  [6]  to  14  SD.  In  this  case,  we  move 
away  from  the  observation  point  of  Case  1  along  the  normal  to  the  triangle  and  compute  the 
integral  at  three  observation  points.  The  integrand  now  is  not  singular  as  in  the  first  four  cases 
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but  the  proximity  of  the  observation  point  to  the  integration  triangle  makes  it  behave  as  if  it 
were.  In  table  4.4,  we  plot  the  difference  between  our  results  and  those  of  Khayat  and  Wilton. 
We  get  excellent  agreement  in  general  except  for  the  real  part  of  Case  5  where  we  note  a 
significant  difference.  We  placed  a  fine  grid  around  this  point  and  found  that  the  Khayat  and 
Wilton  result  falls  off  the  curve.  In  an  exchange  of  correspondence,  the  authors  of  [2]  verified 
the  correctness  of  our  result  using  both  the  original  approach  and  a  modified  one  [8]. 

Table  4. 1 :  Comparison  of  present  approach  to  results  in  references  6  and  7.  Observation  point 
is  inside  the  triangle  and  its  coordinates  are  given  with  respect  to  an  origin 
centered  at  the  right- angle  vertex. 


Observat 
ion  point 
location 

(x,  y,  z) 

in  m 

Case  1 

(0.1,  0.1,  0.0) 

Case  2 

(0.2,  0.2,  0.0) 

Case  3 

(0.3,  0.3,  0.0) 

Case  4 

(0.4,  0.4,  0.0) 

Our 
result 
(7  SD) 

1.89857266168680 

-i 

0.309643085617962 

2.246285006785350 

-i 

0.311143518184329 

2.381002978088960 

-i 

0.311826316300869 

2.283869854898920 

-i 

0.311688243312515 

Our 
result 
(15  SD) 

1.89857266176846 
-i  0.309643085636859 

2.246285006965140 
-i  0.311143518212247 

2.381002978727480 
-i  0.311826316345215 

2.283869855108430 
-i  0.311688243332126 

Khayat 

and 

Wilton 

1.898579 

-i 

0.3096492 

2.246288 

-i 

0.3111477 

2.381009 

-i 

0.3118314 

2.283890 

-i 

0.3116960 

Rossi 

and 

Cullen 
Mathem 
atica 
(6  SD) 

1.89857 

-i 

0.309643 

2.24628 

-i 

0.311144 

2.38099 

-i 

0.311826 

2.28386 

-i 

0.311688 

Khayat 

and 

Wilton 
(14  SD) 

1.89857266176845 

-i 

0.30964308563686 

2.24628500696514 

-i 

0.31114351821225 

2.38100297872747 

-i 

0.31182631634521 

2.28386985510842 

-i 

0.31168824333213 
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Table  4.2:  Difference  between  our  results  (15  SD  rounded  to  14)  and  those  of  Khayat 
and  Wilton  (14  SD):  real  and  imaginary  parts. 


The  points  below  are  for  K&W's 
coordinate  system  (centered  at  the 
right  angle).  Our  system  is  centered 
at  the  centroid. 

Re{K&W  14  SD}  -  Re  {our  method, 
DP,  14  Taylor  Terms} 

Im{K&W  14  SD}  -  Im{our  method, 
DP,  14  Taylor  Terms} 

Case  1:  (.1,  .1,  0) 

0.0000000000001 

0.00000000000000 

Case  2:  (.2,  .2,  0) 

0.0000000000001 

0.00000000000000 

Case  3:  (.3,  .3,  0) 

0.0000000000000 

0.00000000000001 

Case  4:  (.4,  .4,  0) 

0.0000000000000 

0.00000000000000 

Table  4.3:  Observation  point  near  integration  triangle. 
Observation  point  coordinates  (0.1,  0.1,  z). 


Case 

5,  z  =  0.0001 

6,  z  =  0.01 

7,  z  =  0.1 

Re  {/n4i} 

1.897944525246840 

1.837558164829700 

1.429705163246540 

Re  {K&W} 

1.89795445129807 

1.83755816482970 

1.42970516324653 

Im{/(M)} 

-0.309643085431937 

-0.309641036420311 

-0.309438204123196 

Im{K&W} 

-0.30964308543194 

-0.30964103642031 

-0.30943820412320 

Table  4.4:  Differences  for  Cases  5-7  between  our  results  and  those  of  Khayat  and  Wilton. 


Case 

Re  {K&W,  14  SD}  -  Re  {our  method, 
15  SD,  14  Taylor  Terms} 

Im{K&W,  14  SD}  -  Re  {our  method, 
15  SD,  14  Taylor  Terms} 

Case  5:  (.1,  .1,0.0001) 

0.00000992605123 

0.00000000000000 

Case  6:  (.1,  .1,0.01) 

0.00000000000000 

0.00000000000000 

Case  7:  (.1,. 1,0.1) 

-0.00000000000001 

0.00000000000000 

This  concludes  our  validation  tests.  We  will  discuss  validation  further  in  the  second  part  of 
this  report. 


5.  CONCLUSIONS 

We  have  introduced  a  new  method  for  computing  the  inner  integral  of  the  impedance  elements 
in  the  Rao-Wilton-Glisson  [J_]  formulation  of  the  method  of  moments  in  electromagnetics 
(Sections  1  and  2).  The  distinguishing  feature  of  this  method  is  that  it  can  compute  the  integral  to 
a  prescribed  precision.  We  know  of  no  other  method  that  can  do  this.  The  method  is  valid  for  all 
observation  points  that  lie  within  a  sphere  with  center  the  triangle’s  centroid  and  radius  of  one 
half  of  a  wavelength  (Section  3).  This  restriction  is  because  the  two  iteration  formulas  we  use  do 
not  converge  everywhere  in  space  but  have  a  common  domain  of  convergence  in  the  interior  of 
this  sphere.  There  are  also  points  outside  this  sphere  where  we  have  convergence;  they  are  of  no 
consequence,  however,  since  we  can  find  faster  ways  of  obtaining  the  same  accuracy  outside  this 
sphere.  This  is  the  subject  of  the  second  part  of  this  report. 
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APPENDIX  A 

CONVERSION  OF  I,  ,  TO  A  LINE  INTEGRAL 

(") 


Consider  the  integral 

J!l)(r')  =  -(*•/  -p')  4)  (r') 


=  (n  +  \)\\_(x  -xI)x  +  {y' -yi)f\ 


-x'\  +(y-y'\  +  h 1  —r' 


-x'Y  +(y-y')  +  h 2 


-dxdy 


(A.l) 


and  use  the  transformation 


£  =  x-x',  ij  =  y-y' 


(A.2) 


to  write 


JS=(n  +  1)J(cA  +  ^)L^ 


+rj2+h2  -r'J 


2+rj2+h2 


dd,dr] 


(A.  3) 


where 


c,=x'-x,,  b,=y'-yr 


(A.4) 


We  define 


f(l) 

Xn) 


(%,ri)  =  {n  +  \)(clx  +  bly) 


i^%2+ri2+h2  -r'j 

yl?+r,2+h2 


(A.5) 


and  rewrite  it  as  a  complex-valued  function 


Yln\(^d)  =  (n  +  \)(cl+ibl) 


2  +rj2  +h2  —r' 
sle+r+h2 


(A.6) 
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We  make  next  the  transformation 


f=f+/<7.  =■  i=^r~  (a.7) 

2  2i 

to  get 

UcC*+h2-r')n 

VS(C,n  =  (n  +  l)(cl+ib,) V  l—^r  ’  .  (A.  8) 

jcc  +h2 


We  need  the  following  integral 


(«")  Af  V^,t)dT  =  (n  +  \)C-^p-\ 


c,  +ib,  rC 


' r  +  h 2  —r' 


■\ICT  +  b2 


-dr 


(c,  +ib,)^£C+h 

J 


2  t 

—  r 


n+i 


We  proceed  to  define  functions / and  g  as 


(A. 9) 


$>(£,)  =  Re{/“(<r,C)}  =  Re 

-c't+b'y?W+*2 -r')“‘  = 


(  i; )( -Jl1  +T  +lr  -  '■') 


(n+i)(f2+>r) 


£!+<r 


c,4  +  b,rj 
P2 


2  -  7  2  f 

+  /z  —  r 


+//>.  )(  Vf ’+'/'+ *' 


«  +  l  "1 


(„+i)(fJ+r) 


b^-cfl 


e+r 


-r'Y  = 


biC-cgi 

P2 


2  .  /  2  f 

+  /z  -  r 


«+i 


where 


(A.  10) 


(A.  11) 


P  = 


(A.  12) 
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For  (A.5)  we  can  now  write 


(& if)  =  V Q (£  7)  +  £  x  Vg((;>  (£  7) . 


r(/) 


(0, 


(A.  13) 


We  re-adjust  (A.  11)  and  (A.  12)  so  as  to  remove  the  singularity  at  p  =  0.  To  this  end,  we 
compute  the  asymptotic  form  of  these  expressions  as  p  — >  0 .  We  let 


£  =  p  cos  cp ,  rj  =  psin(p 
and  substitute  in  (A.  10)  and  (A.  1 1)  to  get 


'(«) 

r(0 


P 

b,  cos  (p-Cj  sin  cp 
P 


2  f 

—  r 


Vp2 +h 


2  —  r 


n+ 1 


(A- 14) 


(A.  15) 
(A.  16) 


Through  a  direct  calculation,  we  find  that 
c,  cos^  +  h,  sin^? 


'(«) 

r(0 


gm,V)  = 


p 

b,  cos  <p-c,  sin  cp 
P 


(\h\-r')"  +O(l),p^>0 
(|/?|-r')"+1  +O(l),p->0 . 


(A.  17) 
(A.  18) 


We  define 


=  \h\-rT  ,  g 


h,  cos(p-c,  sin^/ 
P 


(A.  19) 


These  two  expressions  satisfy  the  Cauchy-Riemann  conditions.  Thus,  we  can  subtract  them 
from  the  corresponding  ones  in  (A.  10)  and  (A.  1 1)  and,  also,  use  (A.l)  to  generate  functions  that 
do  not  have  a  singularity  at  p  =  0 


.0^)=  p2 


(V p 2 +h2  ->*') 


(A. 20) 


out 


£(„)(£>’ 7)=  pl 


^ Ip 2  +h 


2  —r 


n+ 1 


-W-r'T 


(A. 21) 
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From  [4],  we  have  that,  in  general,  if  v  is  of  the  form  (A.  13),  then 
|  v(x,  y)dxdy  =  J  f(x,y)vds  +  J  g(x,y)tds  . 


(A.22) 


0D 


3D 


where  v  is  the  normal  to  a  triangle’s  side  while  t  is  the  tangent  to  it  and  is  positively  oriented 
with  respect  to  the  nonnal  to  the  triangle.  Both  lie  on  the  triangle’s  plane.  From  (A. 3),  (B20),  and 
(A.21),  we  can  then  write 


f  c,^  +  b,d 

Up'+h'-rf'-m-r')** 

ST  P 

\  J  W  \  / 

■J 


bg-cjn 


dT 


P 


(yjp2+h2  ~r')  ~(\h\-r')n+X 


vds 
ids . 


(A.22) 


We  introduce  rectangular  coordinates 

(vj , ij ,nj,  n  =  z,  v.  =  f.  x  h  ,  j  =  1, 2, 3 


(A. 23) 


on  each  side  of  the  triangle,  with  origin  the  endpoint  that  we  encounter  first  in  going  around  the 
triangle  in  the  positive  sense  with  respect  to  the  normal  n  .  Specifically,  for  side  1,  it  is  the  point 
r2 ;  for  side  2  the  point  r3 ;  and  for  side  3  the  point  r, .  Note  that  these  coordinates  bear  the  index 
of  the  side  they  belong  to.  With  this  notation,  (A.22)  becomes 


f-w-rr 


j= 1  Sj  l  P 

3 


>ds 


-m- 

7=1  S: 


P 


—  r 


(t'l 


ds . 


(A. 24) 


For  the  position  vector  to  the  j- th  side,  we  write 
r  =  r/+l  +  stj ,  j  =  1,2,3 

with  the  indices  running  cyclically.  We  also  write 
r'  =  x'x  +  y'y  +  hz  =  p'  +  hz  ,  r,  =  x,x  +  y,y . 


(A. 25) 


(A. 26) 


From  these  and  (A.24),  we  have  that 
£c ,+rjb ,  =(r-p')-(p'-r/) 
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and 


^[-r/c,  =(x-x')(v'-j/)-(v-v')(x'-x/)  =  z{(r-p')x(p'-r/)]  =  [(p'-r/)x2].(r-p') 

(A. 28) 


We  substitute  these  in  (A. 24) 

J(r-p')'(p'-r,) 


7= 1 


P~ 


2  .  7  2  * 

+  h  -r 


n+ 1 


j[(p'~r/)xf]-(r-Pf) 


P 


2  +A2  —  r' 


(M-'T 


c/.v 


(A. 29) 


We  can  further  write  this  as 


y=i 


Vj  [(r  -  P')  •  (P'  -  r/ )]  +  ij  {[(P'  -  r/ ) x  *]  ’  (r  -  P')} 


P 


2  +h 2  -r' 


■(W 


ds 


We  define 


a./+1  =  aJ+ i«,+i  =  Vi  ~  p'  ’  P=r-p'  =  a  /+I  +  sf 


7+1  '  JI7 


(A. 30) 


(A. 31) 


and  write 


*0  [(r  -  P')  ’  (P'  -  )]  +  h  {[(P'  -  )  X  f  ]  •  I1"  -  P')} 

=  ~Vj  [(r  - P')  ’ a/ ]  - {K  xz]-(r- p')}  =  -a,  {l?,  [(r -  p')  •  a, ]  +  ij  {[a,  x  z]  •  (r  - p')}}  (A.32) 

We  resolve  the  vectors  in  (A.32)  along  the  directions  (a,,z  x  a,,z) 

Vj  [(r  -  p')  •  a,  ]  +  ij  [(a,  x  z)  •  (r  -  p')] 

=  {a;(a/-vy)  +  (zxa/)[(zxa/)-v.]}[(r-p,)-az] 

+  {«/  (a,  ■  tj )  +  {z  x  a,  )[(z  x  d, )  •  f,  ]}  [(a,  x  z)  •  (r  -  p')] 

=  «/  {(«/  •  Vj ) [(r  "  P')  •  a, ]  +  (a,  •  ij ) [(a,  x  z)  •  (r  -  p')]} 
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=  a 


+  (zxa/){[(zxa;)-v.][(r-p')-a;]  +  [(z><a/)-f.][(a/xz)-(r-p')]} 

°i  {(«/  •  Vj  )[(r  "  P')  •  at  ]  +  [a,  •  (z  x  v. )]  [(a,  x  z)  •  (r  -  p')]} 

+  (z  x  a,  ){[(z  x  a, )  •  (f,  x  z)]  [(r  -  p')  •  a,  ]  +  [(z  x  a, )  •  f,  ]  [(a,  x  z)  •  (r  -  p')]J 
«/  {(«/  %•)«,  +[^  •(zxa/)](zxa/)}-(r-p') 

-(zXfl;)J(f.  •«/)«/  +  [(zXa/)-f/.](zXa/)J  •(I'-p') 

7  [^••(r-p')]-(fx«/)pA(r-P')]- 


(A. 33) 


Substitution  in  (A. 32)  gives 


^[(r-p'j'lp'-rjj  +  ^j^p'-rjxzj^r-p'^-a,  [v,  •  (r  - p')]  +  (zx  a, )[£,  •  (r  -  p')] 

(A. 34) 

and  this  in  (A. 30) 


t(0  _  _a 

J(»)  “  a/ 


ir (;  ^ 

j= 1  s,  P 


-(|/?|-r')H+1 


+(fxa/)Zj 


v(r-p') 


m  p 


p2  +  h2  -r'j  -(|A|-r') 


ds 


ds  . 


(A. 35) 


Since  the  original  expression  in  (A. 6)  does  not  have  a  component  along  zxa(,  this  term  must  be 
zero  in  (A. 35).  We  can  show  this  by  using  Stokes’  theorem 


11 


'<  < r  -  pM 


i p 


)n+ 1  i  . 

~(\h\  -f'J  ds  =  Jz-Vx^(r-p')F(/?)]cA 


(A. 36) 


with 


Hp)  = 


[sip2  +h2  -(|/7|-r')”+1 


P 


(A. 37) 


But 


V  x  [(r  -  p')F(p)]  =  VF(p) x  (r  -  p')  =  dF^  Vp x  (r  -  p') 


(A. 38) 
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Vp  =  Vlr  -  pi  =  - — — 
P 


(A. 39) 


Thus,  (A. 38)  is  zero  and,  hence,  the  integral  over  the  triangle  in  (A. 36)  is  zero.  We  can  then 
replace  (A. 35)  by 


j= 1 5|  P 


p2+h2-r'  -(jh  -r')  ds . 


(A.40) 


vj\r-p’)  =  vj\rj+x-p'),  j  =  1,2,3 


(A.41) 


and  (A.40)  takes  the  simpler  form 


Jl.)  (rw  -p  )J - - - 


(A.42) 
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APPENDIX  B 

REGIONS  OF  CONVERGENCE  OF  ITERATION  FORMULAS 

We  begin  with  the  iteration  formula  (2.5)  which  we  write  in  the  fonn 

K+1J=AYnJ+XnJ,  n  =  0,1,2,... 

where 


(B.l) 


^  = 


K,j 


K 


n+lj 


Xnj  = 


Kj 


A  = 


~P 


1 2 


1 

-2  r' 


(B.2) 


According  to  ([9],  p.  371),  the  solution  of  this  difference  equation  is 


n— 1 


r..j=JXj+Z,A'~"x.-u’ 


m= 1 


(B.3) 


The  stability  of  this  solution  is  determined  by  the  homogeneous  part  of  (2.5)  which  we  rewrite  as 
Kn+2J+2r'Kn+lJ  +  p'2Knj  =  0 ,  n  =  0,1,2,...  (B.4) 


Seeking  solutions  of  the  form 


KnJ  =<u" 


(B.5) 


we  end  up  with  the  characteristic  equation 
/u 2  +2r'/u  +  p'2  =  0 
whose  roots  are 

p±  =—r'±sJr '2  - p'2  =-r'±\h\. 

We  thus  have  the  two  solutions 


(B-6) 


(B.7) 


K„  ..  = 


n,j+ 


=  P 


K,  .  = 


=  P 


(B.8) 
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According  to  ([10],  Theorem  7. 2. 9. 6),  a  necessary  and  sufficient  condition  for  the  two  solutions 


to  satisfy  the  stability  condition 

J £ 

lim — =  o  (B.9) 

n->cc  n 

is  that  |//±  |  <  1  and,  if  equality  holds,  then  the  root  must  be  simple.  This  implies  that 

r'-\h\<l,  r'  +  \h\<\  (B.10) 

or  that 

r'  +  \h\<l.  (B.  11) 

For  the  two  roots  to  be  equal,  we  must  have  h  =  0  ,  in  which  case 

H±  =  ~r' ,  h  =  0  (B.12) 

and  for  a  stable  solution  we  must  have 

r'<  1.  (B.13) 

From  this  and  (B.l  1),  we  have  the  two  conditions 

1  —  n'2 

\h\  <  — ,  if  0  <  p’  <  1  and  \h\  >  0 

11  2  11  (B.14) 

0<p'<l  ,  if  \h\  =  0 

We  graph  this  in  figure  B.l. 
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p  prime 


Figure  B.l:  Relation  between  \h\  and  p'  for  stability.  The  stable  region  is  the  one  below 
the  curve  and  bounded  by  the  two  axes.  The  point  (1,  0)  is  excluded. 


We  turn  next  to  the  stability  of  the  iteration  formulas  (2.20),  (2.26)  and  (2.32).  Our  discussion 
is  motivated  by  comments  in  (11 11,  p.  142).  The  homogeneous  part  of  all  three  equations  under 
consideration  is 


rn  +  2  2 

v  ,,  .  = - B  v  . 

m+2’J  m  +  3  7  m’J 


m  =  0,1,2,... . 


(B.l  5) 


This  is  not  an  easy  equation  to  solve  since  the  coefficients  depend  on  the  independent  variable. 
We  can,  however,  proceed  as  follows  in  our  test  for  stability.  We  let 

r  =  ^,  0  <  7  <  1  (B.16) 

m  +  3 

and  treat  /as  a  constant.  We  proceed  as  above  to  obtain  a  characteristic  equation 

p2=rB,2  (B.17) 


with  the  obvious  solutions 

/4  =  ±4yBj  •  (B.18) 
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Since  y<  1,  for  the  solutions  j  and  to  be  stable,  we  need  ffito  be  less  or 

equal  to  one 

0  <  Bj  <  1 .  (B.19) 

We  recall  that  we  obtain  B j  by  dropping  the  nonnal  from  the  observation  point  to  the  line  that 
contains  the y'-th  side  of  the  triangle.  This  distance  is  B  j . 

We  wish  to  point  out  that  (B.19)  does  not  imply  (B.14)  and  conversely.  This  is  unfortunate 
because  both  inequalities  have  to  be  satisfied  simultaneously.  As  examples,  think  of  an 
observation  point  on  the  plane  of  the  triangle  that  lies  on  the  extension  of  one  of  the  sides  of  the 
triangle.  Then,  for  that  side,  B .  =  0  while  r'  can  be  greater  than  one.  If  we  now  consider  an 

observation  point  right  above  the  centroid  of  the  triangle  and  at  a  distance  greater  than  0.5,  then 
(B.14)  is  not  satisfied  while  (B.19)  possibly  could. 
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APPENDIX  C 

AN  ALTERNATE  WAY  OF  COMPUTING  K0  j 

We  assume  the  inequalities  in  (2.49)  and  let 

r  =  B/tana.  (C.l) 

We  then  get 


= 


-u. 


l  [  h  *by+ 2 

sec2  ada 

B 

sec  a  +  \h\ 

1  BJ  J 

J 

1  1 

lfi"by+2 

da 

i  B,  , 
("?) 

cos  <2 1 

{Bj  +  |/?  cosa) 

We  now  let 


u  =  tan 

f  a^ 

,  du  =  — sec^ 

f  a  ' 

l~2) 

2  1 

,  1  +  n2  ,  1  -u2 

da  = - da,  cos  a  = - r- 

2  1  +  w2 


(C.2) 


(C.3) 


and  substitute  in  (C.2) 


tan 


—  tan 


t ,  b , 


K0  =  2 B.  f  L2  ' 

2  B 


BJ 

1|  ‘J±M 


du 


1  +  U  j 


|(i+w2) 


1  +  u 


=  2B  [  y 

1  J  Ji  -4 


tan  —tan 

2  (5, 


(l  +  w2 

\du 

(i-"2)[ 

( BJ  - 

h 

|  w~  +  B  j  +  h 

tan 

f 

1  -l 

—  tan 

2 

tan 

f 

1  -l 

—  tan 

2 

fhb/+2) 

l  bj  J 

r  i  i 

i 

du  2  II 

[  BJ  ) 

J 

tan 

-tan  1 
2 

bj  JJ 

i 

1  —  u  1 +  u _ 

ll  tt  1 

B,-\  / 

|J 

tan 

1 

—  tan 

2 

Vb/+,) 

bj  JJ 

du 


2  Bj  + 
U  H - - - 

B  - 

J 


~  -In  1-w 


itan-1  m+2 

2  B 


1 

—  tan 
2  I  B, 


+  ln  1  +  m 


itan"1  ^+2 
2  B 


—  tan 
2  I  B, 
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The  last  line  in  (C.4)  can  be  written  in  a  variety  of  ways  by  applying  trigonometric  and  inverse 
trigonometric  identities  to  the  numerator.  Finally,  we  remind  that  this  fonnula  is  valid  when  the 
conditions  in  (2.49)  are  satisfied.  We  would  also  have  to  consider  the  other  two  cases. 
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PART  2.  NUMERICAL  EVALUATION  OF  INTEGRALS 
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1.  STATEMENT  OF  PROBLEM  AND  APPROACH 

In  the  first  part  of  this  report,  we  developed  a  method  for  computing  the  integral 

i(0(r')=f_hJ - dS,  f(/)(r)  =  r-r; ,  /  =  1,2,3.  (1.1) 

'  K 


Here 


R  =  yjr'2  -2r'  r  +  r2  =^(x-x')  +(v-v')  +  h1  (1.2) 

with 

r  =  xx  +  yy,  r '  =  x'x  +  y'y  +  hz  (1.3) 

and  the  region  of  integration  T  being  the  triangle  in  figure  1.1.  The  origin  of  coordinates  is  the 
centroid  of  the  triangle,  and  the  latter  lies  on  the  xv-plane.  Boldface  letters  denote  vectors.  The 
same  letters  in  italics  denote  the  magnitudes  of  these  vectors  while  a  caret  over  a  letter  denotes  a 
unit  vector.  The  vector  r,  is  the  position  vector  to  the  /-th  vertex  of  the  triangle,  as  shown  in 
figure  1. 
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In  Part  1,  we  found  it  convenient  to  measure  distances  in  wavelengths  and,  hence,  we  rewrite 
the  original  integral  in  the  form 

I(,V)  =  /t2f — — - dS,  f(/)(r)  =  r -  r,  ,  /  =  1,2,3  (1.4) 

*.  R 

with  all  distances  in  wavelengths  (/t).  We  also  found  it  convenient  to  introduce  the  phase  factor 
exp(-/'2/z7-')  and  write 

I(/)(r')  =  A2e~i2'rr  [ 1  (  )  - dS,  f<Z)(r)  =  r -  r;  ,  /  =  1,2,3.  (1.5) 

•  R 

We  reached  the  conclusion  that  we  can  compute  the  integral  in  (1.5)  to  a  prescribed  precision  by 
approximating  the  exponential  with  its  Maclaurin  series  expansion.  We  can  do  this  for 
observation  points  in  the  interior  of  a  sphere  with  center  the  triangle’s  centroid  and  radius  of  half 
a  wavelength.  We  also  found  that  the  number  of  terms  depends  only  on  the  number  of  significant 
digits  required  and  the  longest  side  of  the  triangle  but  not  on  the  position  of  the  observation  point 
within  the  sphere  or  the  shape  of  the  triangle. 

Since  the  observation  point  may  lie  anywhere  in  space  (but  on  the  three  sides  of  the  triangle), 
we  must  develop  a  method  that  provides  answers  to  (1.5)  for  observation  points  outside  the  half¬ 
wavelength  sphere  and  for  prescribed  accuracy.  We  describe  below  a  procedure  we  developed 
for  doing  show.  Where  we  do  not  use  the  method  of  Part  1,  we  use  cubatures"  to  evaluate  (1.5)  to 
a  prescribed  number  of  significant  digits.  The  general  strategy  is  as  follows. 

1 .  Collect  available  cubatures  for  triangles. 

2.  Develop  a  set  of  testing  triangles. 

3.  For  a  prescribed  accuracy,  develop  a  program  that  provides  timing  information  for  the 
cubatures  and  our  method  in  Part  1  to  run  over  the  set  of  triangles. 

4.  For  observation  points  inside  the  half-wavelength  sphere  but  not  on  the  triangle,  run  our 
method  and  those  cubatures  that  are  faster  than  our  method.  Determine  where  those 
cubatures  provide  the  required  accuracy  by  comparing  their  answer  to  that  of  our  method. 
Based  on  this,  define  a  new  sphere,  concentric  with  the  old  one  but  of  radius  less  or  equal 
to  that  of  the  old  one. 

5.  For  observation  points  outside  the  new  sphere,  order  all  cubatures  according  to  increasing 
size  (number  of  points  employed  in  cubature).  Define  a  convergence  criterion  for  the 
sequence  of  cubatures. 

6.  For  a  judiciously  selected  set  of  observation  points  and  all  triangles,  run  the  sequence  of 
cubatures.  Using  the  results,  define  a  set  of  spheres  concentric  with  the  sphere  in  4, 


2  By  cubature  we  mean  a  numerical  method  for  evaluating  a  surface  integral;  we  reserve  the  word  quadrature  for  a 
numerical  method  for  evaluating  a  line  integral. 
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above.  In  the  space  between  two  successive  spheres,  choose  the  smallest  cubature  that 
has  converged.  This  is  the  cubature  to  be  used  in  this  region  for  the  prescribed  accuracy. 

In  the  following  sections,  we  present  the  details  of  this  plan. 

2.  AVAILABLE  CUBATURES  AND  THEIR  TIMING 

The  cubatures  for  triangles  that  we  use  come  from  three  sources  [J_],  [2],  [3].  They  range  in 
size  from  3  points  to  175  points.  We  have  employed  all  the  cubatures  in  table  5.1  of  [J_]  as  well 
as  a  46-point  one  that  Dr.  Taylor  sent  us.  The  latter  integrates  a  15-deg  polynomial  exactly.  We 
have  also  used  a  175-point  cubature  that  Prof.  H.  Xiao  kindly  supplied  to  us.  It  integrates  a  30- 
deg  polynomial.  The  only  cubature  coming  from  [3]  is  the  7-point  one,  with  a  misprint  corrected. 
We  used  it  because  it  is  a  popular  one  in  the  electromagnetics  community.  Our  guess  is  that  it 
can  at  best  integrate  a  4-deg  polynomial  exactly.  We  exhibit  all  these  cubatures  in  table  2.1  and 
plot  the  relationship  between  polynomial  degree  and  cubature  points  in  figure  2.1. 

Table  2. 1 :  Size  of  cubatures  employed  and  maximum  degree  of  polynomial 

they  integrate  exactly. 


Number  of 
cubature  points 

Maximum  degree  of  polynomial  that  can 
be  integrated  exactly 

3 

2 

6 

4 

7 

4  (?) 

10 

5  i 

15 

7 

21 

9  ! 

28 

ii  ! 

36 

13  1 

45 

14 

46 

15  | 

55 

16 

66 

18  [ 

78 

20  ! 

91 

21  ! 

105 

23 

120 

25  i 

175 

30  ! 
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The  standard  triangle  over  which  these  cubatures  are  defined  is  an  equilateral  triangle  of  side 
one.  Originally,  most  of  them  were  defined  over  a  right  triangle.  We  made  the  appropriate 
transfonnations  to  map  this  triangle  into  our  standard  equilateral  triangle.  We  use  this  triangle  to 
transfonn  a  cubature’s  points  so  that  they  apply  to  any  other,  arbitrarily  shaped  triangle. 

We  demonstrate  our  procedure  using  7  significant  digits  (SD)  as  the  required  accuracy.  The 
number  of  significant  digits  is  a  variable  in  our  programs.  For  demonstration  purposes,  however, 
we  have  chosen  7  SD.  This  number  corresponds  to  single  precision  (32  bits)  in  the  IEEE-754 
floating-point  standard  [4] . 

We  also  set  the  wavelength  A  equal  to  one  meter.  This  does  not  entail  any  loss  of  generality;  it 
simply  allows  for  easier  calculating.  In  our  programs,  the  wavelength  is  a  variable. 

We  require  that  the  longest  side  of  the  triangle  is  no  longer  than  one  tenth  of  a  wavelength.  In 
a  general  MoM  code  that  uses  our  scheme,  the  entire  grid  would  be  scanned  and  the  longest  side 
would  be  found  and  used.  Thus,  the  length  of  the  longest  side  in  our  programs  is  a  variable.  For 
demonstrating  our  approach,  however,  we  set  the  maximum  side  at  less  or  equal  to  a  tenth  of  a 
wavelength. 

Finally,  we  mention  that  all  our  programming  was  done  in  C  programming  language. 
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We  proceed  now  with  the  timing  of  our  method  and  that  of  the  cubatures.  We  only  need  one 
test  triangle  since  the  cubatures  take  the  same  time  over  all  triangles  and  since  the  number  of 
terms  in  our  method  is  fixed  once  the  length  of  the  maximum  side  is  fixed.  The  timing  algorithm 
took  into  account  the  time  to  transfonn  from  the  standard  (equilateral)  triangle  to  the  test  triangle 
in  figure  2.2.  Although  the  two  are  identical  in  this  case,  this  time  must  be  accounted  for  in 
general.  Subroutines  were  written  in  a  way  so  as  to  make  the  timing  comparison  between  our 
method  and  the  cubatures  as  fair  as  possible.  The  results  are  shown  in  figure  2.3.  Time  is 
nonnalized  to  our  method.  All  cubatures  lie  on  a  straight  line  since  time  to  apply  them  is  directly 
proportional  to  a  cubature’s  number  of  points.  We  see  that  cubatures  with  21  or  fewer  points  are 
faster  than  our  method,  while  cubatures  with  28  or  more  points  are  slower.  This  does  not  mean 
that  we  discard  cubatures  that  are  slower.  We  shall  use  them  below  to  establish  convergence  for 
the  sequence  of  cubatures.  Of  immediate  interest,  however,  are  cubatures  that  are  faster  than  our 
method. 


y 


►  x 


(0.5,- 0.5/ V3) 


Figure  2.2:  Equilateral  test  triangle  for  timing  cubatures  and  our  method. 
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Figure  2.3:  Cubature  time  normalized  to  that  of  our  method  (7  SD). 

3.  DIVISION  OF  OBSERVATION  SPACE  INTO  SPHERICAL  REGIONS 

In  the  previous  section  we  determined  that  cubatures  with  2 1  or  fewer  points  are  faster  than 
our  method.  In  this  section  we  resolve  the  question  of  whether  we  can  use  them  within  the  sphere 
of  application  of  our  method  and  obtain  the  required  7-SD  accuracy.  We  also  determine  which 
cubature  to  use  in  the  rest  of  the  observation  space.  To  this  end  we  introduce  a  set  of  25  test 
triangles.  These  triangles  are  shown  in  Appendix  A.  They  have  been  chosen  deterministically 
rather  than  through  a  random  process  that,  say,  generates  triangle  vertices  randomly  subject  to 
the  condition  that  the  longest  side  of  a  triangle  cannot  be  larger  than  A/ 10 .  We  feel  that  our  set 
of  triangles  is  representative  of  those  we  encounter  in  practice,  containing  both  well  behaved  and 
nearly  pathological  ones. 

We  also  generated  a  set  of  approximately  10,000  observation  points  spread  on  a  hemisphere 
above  each  of  these  triangles  and  with  center  the  centroid  of  the  triangle.  This  number  is  the 
same  for  all  25  triangles  once  the  radius  of  the  hemisphere  has  been  fixed.  It  will,  however,  vary 
from  hemisphere  to  hemisphere.  Observation  points  on  a  hemisphere  below  a  triangle  are  not 
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necessary  since  we  have  symmetry  with  respect  to  the  triangle’s  plane.  The  points  are  evenly 
spread  with  respect  to  the  observation  angles.  We  also  used  a  set  of  randomly  chosen  points 
with,  essentially,  the  same  outcome.  In  the  tables  below  we  present  the  results  for  the  set  of 
evenly  spaced  points. 

Using  the  set  of  triangles  and  the  set  of  observation  points,  we  computed  the  three  integrals  in 
(1.5).  These  are  vectors  that  lie  in  the  triangle’s  plane  and  their  components  are  complex 
numbers;  thus,  we  have  3  vectors  of  2  components  each,  with  the  components  a  pair  of  real 
numbers  each,  a  total  of  12  real  numbers.  We  also  computed  the  scalar  integral  that  results  when 
we  eliminate  the  vector  function  from  the  integrand  of  (1.5).  This  is  a  complex  number  and, 
hence,  it  contributes  two  more  real  numbers  for  a  grand  total  of  14. 

We  next  define  what  we  mean  by  average  failure  rate  (AFR)  for  a  specific  cubature.  Suppose 
we  have  the  correct  7-SD  answer  for  every  triangle  and  every  observation  point  that  lies  on  a 
given  hemisphere.  For  one  of  the  triangles  (and  for  the  specific  cubature)  we  compute  the  14  real 
numbers  for  all  10,000  or  so  observation  points.  We  compare  each  such  number  to  the  correct  7- 
SD  answer  by  forming  the  relative  error  and  testing  it  to  determine  whether  it  satisfies  the 
inequality. 


(computed  number)  -  (correct  to  7-SD  number) 
correct  to  7-SD  number 


<  5  10s. 


(3.1) 


We  count  a  failure  if  this  inequality  is  not  satisfied.  The  total  number  of  failures  over  all 
observation  points  (those  that  are  at  the  same  distance  from  the  centroid  of  a  triangle),  and  for  a 
specific  triangle  and  cubature  is  the  failure  rate  (FR).  We  express  this  in  percent.  As  an  example, 
if  we  have  1,500  failures  out  of  10,000  observation  points  tested,  then  the  FR  is 
(1, 500/10, 000)xl00  =  15%. 

For  a  fixed  cubature  and  hemisphere,  we  find  the  FR  for  each  of  the  25  test  triangles.  We  then 
average  the  25  numbers  to  obtain  an  AFR  for  the  specific  hemisphere  and  cubature.  Besides  the 
AFR,  we  also  keep  track  of  the  minimum  FR  (MinFR)  and  maximum  FR  (MaxFR).  These  two 
numbers  come  out  of  two  triangles,  the  one  that  has  the  minimum  FR  and  the  one  that  has  the 
maximum  FR.  We  present  a  sample  of  such  an  output  in  table  3.1.  This  table  is  for  the  7-point 
cubature  (C7)  and  for  all  the  observation  points  on  a  hemisphere  of  radius  r'  =  0.32  .  We  present 
AFR,  MinFR  and  MaxFR  for  the  scalar  integral’s  real  and  imaginary  part.  We  also  present 
results  for  one  of  the  vector  integrals  in  (1.5).  In  this  case,  we  have  four  entries  since  the  vector 
integral  has  two  components,  each  a  complex  number. 
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Table  3.1:  Failure  rates  for  the  7-point  cubature  and  for  the  observation  points  on  a  hemisphere 
of  radius  r'  =  0.32, .  (Isc  real:  real  part  of  scalar  integral;  VI  x  re:  real  part  of  the  x-component  of 

the  vector  integral  in  ( 1 .5)  with  /  =  1 ) 


i  FR 

Isc  real 

Isc  Imag 

VI  x  re 

VI  x  im 

VI  y  re 

VI  y  im 

AFR 

33.04 

0.00 

89.85 

23.90 

93.15 

3.65 

MinFR 

0.00 

0.00 

62.51 

0.00 

82.32 

0.00 

MaxFR 

87.32 

0.00 

99.47 

98.39 

97.33 

36.60 

We  next  describe  how  we  classify  error.  We  are  dealing  with  14  entries  and  there  are  many 
ways  to  define  error.  We  can  work  with  each  of  the  14  entries  separately.  We  can  also  work  with 
all  14  entries  at  once  and  call  a  failure  if  a  specified  number  of  entries  fail  to  produce  7  SD.  Of 
the  many  possibilities,  we  choose  what  we  consider  a  natural  way.  The  scalar  integral  consists  of 
two  entries  (real  and  imaginary);  each  vector  integral  consists  of  four  entries,  since  it  has  two 
complex -valued  components.  We  declare  a  failure  for  the  scalar  integral  if  any  of  its  two  real 
numbers  fails  the  test  (3.1);  similarly,  we  declare  a  failure  for  any  of  the  three  vector  integrals  if 
any  of  its  four  real  numbers  fails  the  test  (3.1).  Under  this  scheme,  we  use  the  largest  number  of 
failures  within  each  of  the  four  groups  of  numbers  to  be  the  FR  of  that  group.  As  an  example, 
table  3.1  would  reduce  to  table  3.2. 

Table  3.2:  Failure  rates  for  the  scalar  and  first  vector  integral  in  table  3.1 
according  to  the  rules  described  above. 


FR 

Isc 

VI 

AFR 

33.04 

93.15 

MinFR 

0.00 

82.32 

MaxFR 

87.32 

99.47 

The  last  issue  we  need  discuss  before  presenting  numerical  results  is  what  we  mean  by 
convergence  of  the  sequence  of  cubatures.  We  look  at  a  single  triangle  and  a  single  observation 
point  and  compute  the  set  of  14  real  numbers  using  all  the  cubatures,  ordered  from  the  smallest 
to  the  largest  size.  For  each  of  the  14  entries  we  compare  how  this  number  changes  as  a  function 
of  cubature  size.  Specifically,  we  compare  successive  results  and  test  them  whether  they  agree  to 
7  SD.  The  test  we  use  is 


C(n,m)-C(n  +  \,m) 


C(n  +  l,m) 


<5-10  8  ,  n  =  1,2,. ..,17; ,  m  =  l,2,...,14 


(3.2) 


where  C(njn)  stands  for  the  result  of  the  nth  cubature  for  the  mth  real  number.  Here,  we  have 

placed  the  17  cubatures  of  table  2.1  in  an  one-to-one  correspondence  with  the  natural  numbers. 
For  each  m,  we  look  for  rt ,  the  smallest  n  that  satisfies(3.2).  If  (3.2)  is  satisfied  for  n  +1  and 
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n  +  2  also,  we  declare  the  sequence  to  have  converged  and  its  value  to  be  that  of  n  +  3  .  In  our 
program,  the  number  of  successive  comparisons  is  a  variable',  thus,  it  can  be  smaller  or  larger 
than  three. 

We  now  present  some  tables  with  results  we  obtained  using  the  triangles,  observation  points, 
cubatures  and  rules  we  described  above.  The  raw  data  sets  we  used  are  too  large  to  present  here. 
As  an  example,  in  Appendix  B  we  present  raw  data  for  the  10-point  cubature  for  a  limited 
number  of  r'  values.  It  is  from  this  kind  of  tables  that  we  obtain  the  tables  that  follow. 

In  table  3.3,  we  present  AFR  results  for  the  scalar  integral  and  for  the  four  cubatures  that  are 
faster  than  our  method.  The  number  of  SD  required  is  7,  as  it  is  throughout  this  discussion.  We 
see  that  cubature  C7  and  CIO  are  not  reliable  at  any  distance  from  the  centroid;  on  the  contrary, 
Cubature  C15  can  be  applied  for  r'>0.2  and  C21  for  r'>0.15  wavelengths.  We  also  see  that 
the  error  starts  increasing  at  large  distances.  We  will  return  to  this  point.  If  we  ignore  the 
behavior  at  large  distances,  then  C15  and  C21  provide  7  SD  with  an  AFR  of  less  than  1%  from 
the  stated  distances  on.  The  contents  of  this  table  appear  in  graphical  fonn  in  figure  3.1. 

Table  3.3:  Scalar  integral  AFR  for  four  cubatures  that  are  faster  than  our  method.  The  number 
we  present  in  each  cell  is  the  larger  of  the  AFR  for  the  real  and  imaginary  parts  of  the  integral. 


r'  (in 

wavelengths) 

IscAFR 

C7 

Isc  AFR 

CIO 

Isc  AFR 

C15 

Isc  AFR 

C21 

!  0.1 

99.61 

99.22 

77.37 

9.85 

0.15 

97.83 

96.74 

7.63 

0.00 

0.2 

91.47 

87.18 

0.00 

0.00 

0.3 

33.05 

21.43 

0.00 

0.00 

0.4 

0.18 

0.00 

0.00 

0.00 

0.5 

97.84 

95.72 

0.01 

0.00 

0.6 

0.00 

0.00 

0.00 

0.00 

0.7 

0.00 

0.00 

0.00 

0.00 

0.8 

0.00 

0.00 

0.00 

0.00 

i  0.9 

0.00 

0.00 

0.00 

0.00 

1 

81.02 

76.96 

0.08 

0.00 

1.5 

77.20 

73.24 

0.06 

0.00 

2 

74.20 

70.70 

0.07 

0.00 

2.5 

72.08 

69.10 

0.07 

0.00 

!  3 

71.15 

67.96 

0.06 

0.00 

5 

71.56 

68.92 

0.05 

0.01 

10 

69.97 

66.90 

0.09 

0.02 

i  25 

70.18 

66.98 

0.08 

0.02 

67.76 

65.67 

0.09 

0.04 

100 

64.87 

64.07 

0.17 

0.11 
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500 

56.83 

59.32 

1.69 

1.36 

1000 

56.92 

60.67 

3.77 

2.99 

Distance  of  Observation  Point  from  Centroid  in  Wavelengths 

Figure  3.1:  Data  of  table  3.3  in  graphical  form:  AFR  for  scalar  integral  versus  distance  of 

observation  point  from  centroid. 

In  Appendix  C.  we  provide  the  same  information  for  the  three  vector  integrals.  Based  on  these 
tables  and  the  criterion  that  the  AFR  must  be  less  than  1%,  we  can  choose  for  each  integral  the 
sphere  with  radius  r'  outside  of  which  we  apply  a  cubature  of  a  certain  size.  The  most  efficient 
(fastest)  way  of  computing  the  four  integrals,  however,  is  a  simultaneous  computation  of  all  four 
(see  Appendix  D);  thus,  it  pays  to  use  the  same  cubature  and  same  radius  r'  for  all  four 
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integrals.  To  this  end,  we  construct  tables  3.4  and  3.5  where,  using  table  3.3  and  the  tables  in 
Appendix  C,  we  choose  the  maximum  AFR  for  each  of  the  cubatures  C15  and  C21  (we  have 
eliminated  the  other  two  as  not  providing  enough  accuracy),  respectively.  We  also  present  these 
results  in  graphical  form  in  figure  3.2.  From  table  3.4,  we  see  that  we  can  use  C15  outside  a 
sphere  of  radius  r'  =  0.4  with  an  AFR  of  less  than  1%  while,  from  table  3.5,  this  radius  becomes 
equal  to  0.2.  In  figure  3.3,  we  display  these  findings  and  we  call  this  conclusion  a  possible 
strategy  for  dividing  the  observation  space. 

Table  3.4:  AFR  for  the  four  integrals  using  C15.  The  last  column  exhibits  the  maximum 

of  the  four  columns  preceding  it. 


r'  (in 

wavelengths) 

Isc  AFR 

C15 

VI  AFR 

C15 

V2  AFR 

C15 

V3  AFR 

Cl  5 

Maximum  AFR 

0.1 

77.37 

94.28 

94.24 

92.41 

94.28 

0.15 

7.63 

47.31 

50.69 

27.08 

50.69  : 

0.2 

0.00 

23.41 

22.30 

6.37 

23.41  : 

0.3 

0.00 

2.62 

2.53 

0.00 

2.62 

0.4 

0.00 

0.18 

0.29 

0.00 

0.29 

0.5 

0.01 

0.19 

0.06 

0.05 

0.19  . 

0.6 

0.00 

0.05 

0.04 

0.00 

0.05 

0.7 

0.00 

0.05 

0.03 

0.00 

0.05  | 

0.8 

0.00 

0.03 

0.02 

0.00 

0.03 

0.9 

0.00 

0.02 

0.02 

0.00 

0.02 

1 

0.08 

0.05 

0.03 

0.03 

0.08 

1.5 

0.06 

0.07 

0.01 

0.01 

o.o7  ; 

2 

0.07 

0.00 

0.00 

0.00 

0.07  j 

2.5 

0.07 

0.00 

0.00 

0.00 

0.07  ; 

3 

0.06 

0.00 

0.00 

0.00 

0.06 

5 

0.05 

0.00 

0.00 

0.00 

0.05  | 

10 

0.09 

0.00 

0.00 

0.00 

0.09  ! 

25 

0.08 

0.00 

0.00 

0.00 

0.08  : 

50 

0.09 

0.00 

0.00 

0.00 

0.09  ; 

100 

0.17 

0.00 

0.00 

0.00 

0.17 

500 

1.69 

0.00 

0.00 

0.00 

1.69 

1000 

3.77 

0.00 

0.00 

0.00 

3.77 
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Table  3.5:  AFR  for  the  four  integrals  using  C21.  The  last  column  exhibits  the  maximum 

of  the  four  columns  preceding  it. 


r'  (in 

wavelengths) 

Isc  AFR 

C21 

VI  AFR 

C21 

V2  AFR 

C21 

V3  AFR 

C21 

Maximum  AFR 

0.1 

9.85 

55.88 

48.52 

30.45 

55.88  ! 

0.15 

0.00 

4.26 

5.32 

0.01 

5.32  j 

0.2 

0.00 

0.34 

0.57 

0.00 

0.57 

CO 

o 

0.00 

0.02 

0.02 

0.00 

0.02  ; 

0.4 

0.00 

0.00 

0.00 

0.00 

0.00  ' 

0.5 

0.00 

0.00 

0.00 

0.00 

o.oo  ; 

0.6 

0.00 

0.00 

0.00 

0.00 

0.00  | 

0.7 

0.00 

0.00 

0.00 

0.00 

0.00 

CO 

o 

0.00 

0.00 

0.00 

0.00 

0.00 

0.9 

0.00 

0.00 

0.00 

0.00 

0.00 

1 

0.00 

0.00 

0.00 

0.00 

0.00  ; 

1.5 

0.00 

0.00 

0.00 

0.00 

0.00  | 

2 

0.00 

0.00 

0.00 

0.00 

0.00 

2.5 

0.00 

0.00 

0.00 

0.00 

0.00 

3 

0.00 

0.00 

0.00 

0.00 

0.00 

5 

0.01 

0.00 

0.00 

0.00 

0.01  i 

10 

0.02 

0.00 

0.00 

0.00 

0.02 

25 

0.02 

0.00 

0.00 

0.00 

0.02  ; 

0.04 

0.00 

0.00 

0.00 

0.04 

100 

0.11 

0.00 

0.00 

0.00 

0.11 

500 

1.36 

0.00 

0.00 

0.00 

1.36  i 

1000 

2.99 

0.00 

0.00 

0.00 

2.99  | 
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Distance  of  Observation  Point  from  Centroid  in  Wavelengths 

Figure  3.2:  Graphical  form  of  right-hand  columns  of  tables  3.4  and  3.5:  Maximum  AFR  as  a 
function  of  distance  of  observation  point  from  centroid. 
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Figure  3.3:  A  possible  zoning  strategy  for  7  SD  based  on  a  less  than  1%  AFR. 

Another  possible  strategy  for  dividing  the  observation  space  is  to  use  the  MaxFR.  In  Appendix 
E,  we  present  the  relevant  tables  and  use  them  to  construct  tables  3.6  and  3.7,  and  we  also 
display  their  contents  graphically  in  figure  3.4.  Using  these  tables,  we  present  in  figure  3.5 
another  possible  strategy  based  on  a  MaxFR  of  less  than  1%. 
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Table  3.6:  MaxFR  for  the  four  integrals  using  Cl 5.  The  last  column  exhibits  the  maximum 

of  the  four  columns  preceding  it. 


r'  (in 

wavelengths) 

Isc  MaxFR 

Cl  5 

VI  MaxFR 
C15 

V2  MaxFR 

Cl  5 

V3  MaxFR 

Cl  5 

Maximum 

MaxFR 

0.1 

96.40 

99.79 

99.62 

97.53 

99.79  ! 

0.15 

26.46 

97.46 

97.46 

79.99 

97.46 

0.2 

0.00 

76.93 

81.28 

29.22 

81.28 

CO 

o 

0.00 

15.12 

15.92 

0.00 

15.92  ; 

0.4 

0.00 

2.73 

7.19 

0.00 

7.19 

0.5 

0.03 

2.78 

1.45 

0.14 

2.78  | 

0.6 

0.00 

0.59 

0.67 

0.00 

0.7 

0.00 

1.15 

0.25 

0.00 

1.15 

CO 

o 

0.00 

0.28 

0.34 

0.00 

0.9 

0.00 

0.16 

0.45 

0.00 

1 

0.22 

1.18 

0.11 

0.06 

1.18 

1.5 

0.18 

1.82 

0.03 

0.02 

1.82  | 

2 

0.17 

0.03 

0.03 

0.02 

0.17 

2.5 

0.17 

0.01 

0.01 

0.01 

3 

0.15 

0.02 

0.02 

0.02 

5 

0.26 

0.04 

0.07 

0.02 

10 

0.31 

0.04 

0.09 

0.02 

0.31 

25 

0.14 

0.02 

0.04 

0.01 

0.19 

0.01 

0.08 

0.01 

100 

0.34 

0.03 

0.09 

0.04 

500 

2.80 

0.01 

0.04 

0.01 

1000 

6.26 

0.04 

0.09 

0.04 

6.26 
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Table  3.7:  MaxFR  for  the  four  integrals  using  C21.  The  last  column  exhibits  the  maximum 

of  the  four  columns  preceding  it. 


r'  (in 

wavelengths) 

Isc  MaxFR 
C21 

VI  MaxFR 
C21 

V2  MaxFR 
C21 

V3  MaxFR 
C21 

Maximum 

MaxFR 

0.1 

54.96 

97.44 

94.77 

65.99 

97.44  ! 

0.15 

0.00 

18.31 

31.32 

0.26 

31.32  | 

0.2 

0.00 

2.06 

5.03 

0.00 

CO 

o 

0.00 

0.11 

0.30 

0.00 

0.4 

0.00 

0.05 

0.01 

0.00 

0.5 

0.01 

0.01 

0.02 

0.00 

0.6 

0.00 

0.01 

0.00 

0.00 

0.01  ; 

0.7 

0.00 

0.00 

0.00 

0.00 

0.00 

CO 

o 

0.00 

0.00 

0.00 

0.00 

0.00 

0.9 

0.00 

0.00 

0.00 

0.00 

0.00  : 

1 

0.01 

0.00 

0.00 

0.01 

o.oi  ; 

1.5 

0.01 

0.01 

0.00 

0.00 

0.01  | 

2 

0.01 

0.00 

0.00 

0.00 

0.01 

2.5 

0.02 

0.00 

0.00 

0.00 

0.02 

3 

0.02 

0.00 

0.00 

0.00 

0.02 

5 

0.02 

0.00 

0.00 

0.00 

0.02  | 

10 

0.06 

0.00 

0.00 

0.00 

0.06 

25 

0.06 

0.00 

0.00 

0.00 

0.06  ; 

0.13 

0.00 

0.00 

0.01 

0.13 

100 

0.25 

0.01 

0.01 

0.00 

0.25 

500 

2.34 

0.01 

0.02 

0.01 

2.34 

1000 

5.62 

0.01 

0.01 

0.01 

5.62  | 
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Distance  of  Observation  Point  from  Centroid  in  Wavelengths 

Figure  3.4:  Graphical  form  of  right-hand  columns  of  tables  3.6  and  3.7:  Maximum  MaxAFR  as 
a  function  of  distance  of  observation  point  from  centroid. 
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Figure  3.5:  A  possible  zoning  strategy  for  7  SD  based  on  a  less  than  1%  MaxFR. 

This  is  in  brief  the  strategy  for  determining  the  spherical  regions  where  we  employ  our  method 
(see  Part  1  of  this  report)  and  the  C21  and  C15  cubatures.  The  criteria  we  used  in  determining 
these  regions  are  not  unique  and  can  be  altered  according  to  the  preferences  of  the  user  and/or 
the  requirements  of  the  problem.  One  issue  we  have  not  resolved  is  the  rise  in  the  FR  for  the 
scalar  integral  as  the  observation  point  tends  to  infinity.  This  occurs  only  in  the  imaginary  part  of 
the  integral  and  reduces  the  accuracy  of  that  part  to  5  SD,  on  average.  It  is  the  result  of  the 
subtraction  of  two  numbers  that  are  accurate  to  7  SD  and  whose  leading  two  digits  are  the  same; 
thus,  what  remains  after  subtraction  is  accurate  to  5  SD.  This  imaginary  part  is  much  smaller 
than  the  real  part,  by  two  or  more  orders.  We  do  not  know  how  the  loss  of  SD  in  the 
considerably  smaller  part  affects  the  accuracy  of  the  solution  of  the  system  of  equations,  nor  do 
we  think  that  it  is  easy  to  analyze. 


68 


NAWCADPAX/TR-2008/227 


4.  CONCLUSIONS 

For  the  Rao-Wilton-Glisson  formulation  of  the  method  of  moments  [5],  we  have  developed  a 
procedure  for  dividing  the  observation  space  into  spherical  regions  when  computing  the  inner 
surface  integrals.  This  procedure  guarantees  a  prescribed  number  of  SD  in  the  evaluation  of 
these  integrals.  In  the  innermost  region,  we  use  the  method  we  developed  in  Part  1  of  this  report. 
In  the  remaining  two  regions,  we  use  two  different  size  cubatures.  The  radius  of  each  spherical 
region  and  the  cubature  size  depend  on  the  number  of  SD  that  we  require.  The  number  of 
spherical  regions  is  not  fixed  but  can  be  adjusted  by  the  user.  In  our  experience,  three  spherical 
regions  are  sufficient  to  cover  the  entire  observation  space. 

We  can  also  express  the  original  integral  in  terms  of  line  integrals  around  the  triangle’s 
boundary.  These  integrals  can  be  computed  using  a  Gauss-Legendre  quadrature.  Convergence  to 
a  desired  accuracy  can  be  determined  in  the  same  way  as  above  by  using  a  sequence  of  such 
quadratures.  We  have  not  tried  this  approach  but  it  would  be  interesting  to  find  out  how  it 
compares  time  wise  with  the  present  one. 
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APPENDIX  A 

SET  OF  TEST  TRIANGLES 

In  this  appendix  we  present  the  set  of  25  test  triangles  we  use  in  Section  3.  They  appear  in 
figures  A.l  -  A. 4.  In  figure  A.l  we  have  a  right  triangle  whose  hypotenuse  is  equal  to  0.1 
wavelengths.  The  angle  #  starts  at  10°  and  is  given  increments  of  5°,  up  to  45°  for  a  total  of  eight 
triangles. 


Figure  A.l:  Set  of  right  triangles.  AC  =  A/ 10,  0  =  10°  (5°)  45° 

In  figure  A. 2  we  have  a  set  of  obtuse  triangles  whose  maximum  side  is  A  HO .  Again,  the  angle  6 
starts  at  10°  and  is  given  increments  of  5°,  up  to  45°  for  a  total  of  eight  triangles. 
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Figure  A. 2:  Set  of  obtuse  triangles.  AC  =  / 1/10,  AB  =  / 1/20,  9  =  10°  (5°)  45° 

In  figure  A. 3,  we  have  a  set  of  obtuse  isosceles  triangles  with  the  equal  sides  being  XI 20  in 
length.  Again,  the  angle  6  starts  at  10°  and  is  given  increments  of  5°,  up  to  45°  for  a  total  of  eight 
triangles. 


Figure  A. 3:  Set  of  isosceles  triangles.  AB  =  AC  =  XI  20,  9  =  10°  (5°)  45° 
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The  last  of  the  25  triangles  is  an  equilateral  triangle  whose  side  is  equal  to  T/10.  The  triangle 
is  rotated  clockwise  about  its  centroid,  through  an  angle  of  20° . 


Figure  A. 4:  Single  equilateral  triangle  with  side  /L/10 ,  rotated  clockwise  by  20° . 
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APPENDIX  B 

RAW  DATA  FOR  10-POINT  CUBATURE 

Table  B.  1 :  AFR  for  the  10-Point  Cubature 


t 

r 

V1 x re 

V1 x im 

V1 y re 

V1 y im 

1 

BISS 

US 

bus 

Eras 

ffiRISi 

;  0.1 

99.22 

0.00 

99.83 

13.37 

99.83 

0.00 

99.82 

7.33 

99.84 

0.00 

99.78 

0.00 

99.65 

0.15 

96.74 

0.00 

99.19 

16.27 

99.07 

0.00 

99.02 

8.84 

99.08 

0.00 

98.66 

0.00 

98.38 

j  0.2 

87.18 

0.00 

97.97 

18.20 

98.45 

0.00 

98.85 

12.66 

98.51 

0.00 

98.62 

0.00 

96.09 

0.3 

21.43 

0.00 

87.02 

22.97 

89.95 

0.00 

91.20 

21.60 

90.78 

0.00 

89.19 

2.53 

74.76 

;  0.4 

0.00 

0.00 

37.63 

36.29 

30.33 

21.77 

40.57 

44.70 

33.24 

25.33 

21.02 

27.69 

7.55 

i  0.5 

0.00 

95.72 

29.12 

99.39 

7.17 

98.08 

27.39 

99.49 

10.12 

98.10 

7.10 

99.43 

0.00 

Table  B.2:  MinFR  for  the  10-Point  Cubature 


t 

r 

lsc real 

lsc lmag 

V 1 x re 

V1 x im 

V1 y re 

V1 y im 

V2 x re 

V2 x im 

V2 y re 

V2 y im 

V3 x re 

V3 x im 

V3 y re 

V3 y im 

'  0.1 

97.24 

99.38 

99.64 

99.44 

99.64 

99.52 

98.93 

0.00 

91.46 

97.65 

98.24 

96.68 

98.24 

96.72 

96.08 

0.00 

0.2 

63.78 

0.00 

94.23 

0.00 

95.96 

0.00 

96.21 

0.00 

95.96 

0.00 

96.25 

0.00 

91.94 

CO 

o 

0.00 

0.00 

44.78 

0.00 

72.40 

0.00 

66.89 

0.00 

72.40 

0.00 

67.11 

0.00 

40.47 

j  0.4 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

0.00 

’  0.5 

0.00 

92.10 

0.00 

96.48 

0.00 

91.90 

0.00 

97.90 

0.00 

91.90 

0.00 

97.90 

0.00 

Table  B.3:  MaxFR  for  the  10-Point  Cubature 


Iscjmag 

V 1  _x_re 

V1_x_im 

V1_y_re 

V1_y_im 

V2  x  re 

V2_x_im 

V2  y  re 

V2_y_im 

V3_x_re 

V3_x_im 

V3_y_re 

V3_y_im 

0.1 

99.88 

0.00 

100.00 

91.52 

99.97 

0.00 

100.00 

73.55 

99.97 

0.00 

99.95 

0.00 

99.94 

0.00 

0.15 

99.10 

0.00 

99.96 

94.19 

99.74 

0.00 

99.96 

82.07 

99.74 

0.00 

99.71 

0.00 

99.63 

0.2 

96.87 

0.00 

99.83 

95.41 

99.49 

0.00 

99.96 

88.51 

99.51 

0.00 

99.70 

0.00 

99.13 

0.3 

86.51 

0.00 

99.23 

96.65 

96.16 

0.00 

99.41 

95.68 

97.59 

0.00 

97.15 

63.14 

94.11 

0.4 

0.00 

0.00 

99.16 

98.90 

76.28 

66.70 

98.77 

98.68 

87.77 

89.07 

87.80 

80.00 

50.30 

0.5 

0.00 

98.40 

98.58 

99.99 

58.20 

99.54 

98.66 

99.99 

73.73 

99.54 

74.59 

99.99 

0.00 
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APPENDIX  C 

AVERAGE  FAILURE  RATE  FOR  VECTOR  INTEGRALS 

In  tables  C.l  -  C.3,  we  present  AFR  results  for  the  three  vector  integrals  and  for  the  same 
cubatures  as  in  table  3.3. 


Table  C.  1 :  First  vector  integral  AFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  AFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

VI  AFR 

VI  AFR 

VI  AFR 

VI  AFR 

C7 

CIO 

C15 

C21  ' 

0.1 

99.90 

99.83 

94.28 

55.88  ' 

0.15 

99.42 

99.19 

47.31 

4.26 

0.2 

98.96 

98.45 

23.41 

0.34 

0.3 

93.15 

89.95 

2.62 

0.02 

0.4 

42.50 

37.63 

0.18 

0.00 

0.5 

99.62 

99.39 

0.19 

0.00  j 

0.6 

35.66 

33.10 

0.05 

0.00  ' 

0.7 

47.05 

40.89 

0.05 

0.00 

0.8 

40.31 

36.34 

0.03 

0.00  | 

0.9 

35.33 

31.18 

0.02 

0.00  ; 

1 

87.61 

81.66 

0.05 

0.00  ! 

1.5 

78.03 

71.95 

0.07 

o.oo  ! 

2 

70.69 

66.11 

0.00 

o.oo  : 

2.5 

66.64 

60.82 

0.00 

o.oo  ; 

3 

64.00 

58.70 

0.00 

0.00 

5 

57.64 

54.42 

0.00 

0.00 

10 

55.40 

51.62 

0.00 

0.00 

25 

57.28 

53.95 

0.00 

0.00 

50 

55.42 

52.91 

0.00 

0.00  1 

100 

55.33 

51.40 

0.00 

0.00  | 

500 

55.40 

52.92 

0.00 

0.00  ; 

1000 

55.32 

51.42 

0.00 

0.00 
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Table  C.2:  Second  vector  integral  AFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  AFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

V2  AFR 

V2  AFR 

V2  AFR 

V2  AFR 

C7 

CIO 

C15 

C21 

0.1 

99.90 

99.84 

94.24 

48.52 

0.15 

99.41 

99.08 

50.69 

5.32 

99.16 

98.85 

22.30 

0.57 

0.3 

93.97 

91.20 

2.53 

0.02 

0.4 

51.43 

44.70 

0.29 

0.00 

0.5 

99.60 

99.49 

0.06 

0.00  : 

0.6 

38.10 

34.28 

0.04 

0.00 

0.7 

44.82 

39.88 

0.03 

0.00  i 

0.8 

44.43 

38.96 

0.02 

0.00 

0.9 

32.52 

27.29 

0.02 

o.oo  : 

1 

87.57 

81.86 

0.03 

0.00  ! 

1.5 

78.21 

72.17 

0.01 

o.oo  ! 

2 

70.70 

66.19 

0.00 

o.oo  : 

2.5 

66.63 

60.85 

0.00 

0.00 

3 

64.00 

58.70 

0.00 

0.00 

5 

57.64 

54.41 

0.00 

0.00 

10 

55.38 

51.62 

0.00 

o.oo  | 

25 

57.27 

53.94 

0.00 

0.00 

50 

55.42 

52.91 

0.00 

0.00 

100 

55.31 

51.40 

0.00 

0.00  1 

500 

55.40 

52.92 

0.00 

0.00  ! 

1000 

55.32 

51.41 

0.00 

o.oo  ! 
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Table  C.3:  Third  vector  integral  AFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  AFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

V3  AFR 

V3  AFR 

V3  AFR 

V3  AFR 

C7 

CIO 

C15 

C21 

0.1 

99.852 

99.78 

92.41 

30.45 

0.15 

99.094 

98.66 

27.08 

0.01 

0.2 

98.964 

98.62 

6.37 

0.00 

0.3 

92.431 

89.19 

0.00 

0.00 

0.4 

34.798 

27.69 

0.00 

0.00 

0.5 

99.555 

99.43 

0.05 

0.00  : 

0.6 

18.483 

14.41 

0.00 

0.00 

0.7 

26.936 

21.09 

0.00 

0.00 

0.8 

24.980 

20.04 

0.00 

0.00 

0.9 

13.323 

10.15 

0.00 

o.oo  : 

1 

87.495 

81.76 

0.03 

0.00  1 

1.5 

78.102 

72.13 

0.01 

0.00  I 

2 

70.724 

66.18 

0.00 

o.oo  ; 

2.5 

66.609 

60.87 

0.00 

0.00 

3 

64.021 

58.55 

0.00 

0.00 

5 

57.629 

54.39 

0.00 

0.00 

10 

55.386 

51.66 

0.00 

0.00  1 

25 

57.270 

53.93 

0.00 

0.00  ; 

50 

55.451 

52.91 

0.00 

0.00 

100 

55.317 

51.41 

0.00 

0.00  ! 

500 

55.393 

52.90 

0.00 

0.00  ! 

1000 

55.320 

51.41 

0.00 

0.00  ! 
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APPENDIX  D 

COMPUTATION  OF  THE  FOUR  INTEGRALS  USING  CUBATURES 

We  show  how  we  computed  the  integrals  (1.1)  as  well  as  the  scalar  integral 

-ikR 

4  =  f- - ds  . 

sc  J  R 


(D.l) 


This  last  integral  is  computed  via  a  cubature  of  N  points 


4=Z 


N  Q~ikR« 


=  ZjW« 

n= 1 


R. 


(D.2) 


We  compute  the  integrand  at  the  N  points  and  store  the  information.  The  x-  components  of  the 
integrals  in  ( 1 . 1)  are  then  computed  via 


ly  . 

I?)  =Zw„(w  -*/)- 


-ikR„ 


R 


(D.3) 


while  the  v-components  by 


I.i,)  =  2X(t„  ~y>)- 


-ikR. , 


n= 1 


R„ 


(D.4) 


77 


NAWCADPAX/TR-2008/227 


APPENDIX  E 

MAXIMUM  FAILURE  RATE  FOR  THE  FOUR  INTEGRALS 

In  tables  E.  1  -  E.4,  we  present  MaxFR  results  for  the  three  vector  integrals  and  for  the  same 
cubatures  as  in  table  3.3. 


Table  E.  1 :  Scalar  integral  MaxFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  MaxFR  for  the  real  and  imaginary  parts  of  the 

integral. 


r'  (in  wavelengths) 

Isc  MaxFR 

Isc  MaxFR 

Isc  MaxFR 

Isc  MaxFR 

C7 

CIO 

C15 

C21 

'  0.1 

99.92 

99.88 

96.40 

54.96 

;  0.15 

99.31 

99.10 

26.46 

0.00 

0.2 

98.34 

96.87 

0.00 

0.00 

CO 

o 

87.32 

86.51 

0.00 

0.00 

!  0.4 

4.40 

0.00 

0.00 

0.00 

:  o.5 

98.69 

98.40 

0.03 

0.01 

,  0.6 

0.00 

0.00 

0.00 

0.00 

:  o.7 

0.00 

0.00 

0.00 

0.00 

i  0.8 

0.00 

0.00 

0.00 

0.00 

0.9 

0.00 

0.00 

0.00 

0.00 

1 

93.67 

92.49 

0.22 

0.01 

1.5 

92.21 

91.44 

0.18 

0.01 

;  2 

90.35 

89.52 

0.17 

0.01 

:  2.5 

87.45 

85.66 

0.17 

0.02 

3 

90.97 

90.11 

0.15 

0.02 

5 

91.53 

89.52 

0.26 

0.02 

10 

94.14 

93.98 

0.31 

0.06 

:  25 

85.15 

84.98 

0.14 

0.06 

90.82 

89.71 

0.19 

0.13 

!  100 

86.40 

87.97 

0.34 

0.25 

j  500 

84.21 

85.47 

2.80 

2.34 

;  iooo 

83.33 

86.03 

6.26 

5.62 
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Table  E.2:  First  vector  integral  MaxFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  MaxFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

VI  MaxFR 

VI  MaxFR 

VI  MaxFR 

VI  MaxFR 

C7 

CIO 

C15 

C21 

0.1 

99.99 

100.00 

99.79 

97.44 

!  0.15 

99.98 

99.96 

97.46 

18.31 

!  0.2 

99.88 

99.83 

76.93 

2.06 

i  0.3 

99.47 

99.23 

15.12 

0.11 

!  0.4 

99.48 

99.16 

2.73 

0.05 

,  0.5 

100.00 

99.99 

2.78 

0.01 

0.6 

97.43 

96.30 

0.59 

0.01 

!  0.7 

97.98 

97.11 

1.15 

0.00 

0.8 

96.39 

93.06 

0.28 

0.00 

J  0.9 

96.19 

94.58 

0.16 

0.00 

1 

96.25 

94.41 

1.18 

0.00 

1.5 

94.90 

92.85 

1.82 

0.01 

i  2 

88.93 

87.21 

0.03 

0.00 

2.5 

88.10 

85.90 

0.01 

0.00 

3 

86.88 

84.28 

0.02 

0.00 

'  5 

84.73 

82.60 

0.04 

0.00 

10 

87.32 

85.56 

0.04 

0.00 

25 

84.14 

82.00 

0.02 

0.00 

!  50 

84.92 

82.33 

0.01 

0.00 

,  100 

87.50 

85.75 

0.03 

0.01 

!  500 

84.98 

82.38 

0.01 

0.01 

,  1000 

87.50 

85.77 

0.04 

0.01 
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Table  E.3:  Second  vector  integral  MaxFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  MaxFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

V2  MaxFR 

V2  MaxFR 

V2  MaxFR 

V2  MaxFR 

C7 

CIO 

C15 

C21 

0.1 

100.00 

100.00 

99.62 

94.77 

!  0.15 

99.98 

99.96 

97.46 

31.32 

!  0.2 

99.93 

99.96 

81.28 

5.03 

99.50 

99.41 

15.92 

0.30 

0.4 

99.23 

98.77 

0.01 

,  0.5 

100.00 

99.99 

1.45 

0.02 

0.6 

92.61 

91.32 

0.67 

0.00 

!  0.7 

97.18 

96.09 

0.25 

0.00 

0.8 

96.05 

94.59 

0.34 

0.00 

J  0.9 

93.45 

91.25 

0.45 

0.00 

1 

95.87 

94.46 

0.11 

0.00 

1.5 

92.65 

91.83 

0.00 

i  2 

88.90 

87.22 

0.03 

0.00 

2.5 

87.84 

85.91 

0.01 

0.00 

3 

86.56 

84.23 

0.02 

0.00 

'  5 

84.62 

82.36 

0.07 

0.00 

10 

87.45 

85.67 

0.09 

0.00 

25 

84.07 

81.92 

0.04 

0.00 

!  50 

85.01 

82.28 

0.08 

0.00 

,  100 

87.50 

85.67 

0.09 

0.01 

!  500 

85.01 

82.29 

0.04 

0.02 

,  1000 

87.52 

85.71 

0.09 

0.01 
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Table  E.4:  Third  vector  integral  MaxFR  for  four  cubatures  that  are  faster  than  our  method.  The 
number  we  present  in  each  cell  is  the  larger  of  the  MaxFR  for  the  four  real  numbers  that  make  up 
the  real  and  imaginary  parts  of  the  two  components  of  the  integral. 


r'  (in  wavelengths) 

V3  MaxFR 

V3  MaxFR 

V3  MaxFR 

V3  MaxFR 

C7 

CIO 

C15 

C21 

0.1 

99.95 

99.95 

97.53 

65.99 

!  0.15 

99.81 

99.71 

79.99 

0.26 

!  0.2 

99.76 

99.70 

29.22 

0.00 

97.30 

97.15 

0.00 

0.00 

0.4 

88.06 

87.80 

0.00 

0.00 

,  0.5 

99.96 

99.99 

0.14 

0.00 

0.6 

61.53 

57.66 

0.00 

0.00 

!  0.7 

74.29 

68.55 

0.00 

0.00 

0.8 

76.65 

68.62 

0.00 

0.00 

J  0.9 

28.76 

27.16 

0.00 

0.00 

1 

95.60 

94.33 

0.06 

0.01 

1.5 

92.65 

91.83 

0.02 

0.00 

i  2 

88.66 

87.15 

0.02 

0.00 

2.5 

87.72 

85.83 

0.01 

0.00 

3 

86.73 

84.30 

0.02 

0.00 

'  5 

84.64 

82.30 

0.02 

0.00 

10 

87.48 

85.70 

0.02 

0.00 

25 

84.13 

81.87 

0.01 

0.00 

!  50 

84.89 

82.30 

0.01 

0.01 

,  100 

87.42 

85.73 

0.04 

0.00 

!  500 

84.98 

82.31 

0.01 

0.01 

,  1000 

87.42 

85.76 

0.04 

0.01 
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